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Abstract 

These notes summarize a series of works related to the numerical approximation 
of plasma fluid problems. We construct so-called 'Asymptotic-Preserving' schemes 
which are valid for a large range of values (from very small to order unity) of the 
dimensionless parameters that appear in plasma fluid models. Specifically, we are 
interested in two parameters, the scaled Debye length which quantifies how close to 
quasi-neutrality the plasma is, and the scaled cyclotron period, which is inversely 
proportional to the magnetic field strength. We will largely focus on the ideas, in 
order to enable the reader to apply these concepts to other situations. 
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Introduction 



Numerical resolution of perturbation problems 

These notes are about the apphcation of the 'Asymptotic-Preserving' methodology to 
construct schemes for plasma fluids problem which sustain large variations of some of the 
characteristic dimensionless parameters of the plasma. We will be specifically concerned 
with two of these parameters, the scaled Debye length on the one hand and the scaled 
cyclotron period on the other hand. The former quantifies how close to quasi-neutrality 
the plasma is while the latter measures the confinement effects due to the magnetic field. 

Let us consider a singular perturbation problem whose solutions converge to those 
of a limit problem P° when the perturbation parameter e tends to zero. Usually, when 
e <C 1, standard numerical methods (like e.g. explicit methods in the case of time- 
dependent problems) break down. The reason is that the stability condition limits the 
allowed time step to a maximal value which depends on e and tends to zero when e — )■ 0. 
In the case of hyperbolic problems, this problem occurs if one of the wave-speeds tends 
to infinity with e. For instance, in the low Mach-number limit of compressible fluids, the 
acoustic wave speeds tend to infinity when the Mach-number tends to zero. 

In order to overcome such problems, the usual strategy consists in solving the limit 
problem P° instead of P^. For instance, in the case of the Low Mach-number limit, the 
incompressible Euler or Navier-Stokes equations will be solved. However, there are several 
difficulties with this strategy, which we now outline. The first one is that it supposes 
that P° has been previously determined and the second one is that it assumes that 
easy to solve. Both assumptions are by no means obvious. There are cases where the 
determination of the limit problem is difficult if not impossible. Even if P" is well-known, 
it usually involves equations of mixed type (for instance, the small Mach-number flow 
model involves a combination of equations of hyperbolic and elliptic character), where 
constraints (such as the divergence free constraint on the velocity) need to be enforced. 
The abundant literature on the Stokes and Navier-Stokes problem shows that enforcing 
this constraint is not an easy problem. 

The problem complexifies even more when the parameter e is not uniformly small. 
This sentence may sound a little awkward, since e is a number which should have a 
definite fixed value. However, e is usually a ratio of characteristic lengths which may vary 
in space and time. For instance, in a plasma sheath, where the density drops by orders of 
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magnitude, the value of the Debye length changes dramatically. In other cases, such as in 
boundary layers, one must change the scaling length from say the size of the experiment, 
to the typical dimension of the boundary layer. Therefore, the definition of a uniform 
value for e is difficult, and it is more appropriate to view e as a local quantity. 

In such a situation, e may be small in some areas and order unity in other regions. 
Then, the use of the limit problem leads to wrong results in the regions where e is not 
small. To make a more accurate simulation, it is necessary to decompose the simulation 
domain into regions where e = 0(1) and regions where e <^ 1 and to solve in the 
former and P° in the latter. However, the practical realization of this coupling strategy 
is very complex. 

Indeed, first, it requires to define the location where the transition from to P° or 
vice-versa takes place. This is not an obvious question when e has a smooth rather than 
abrupt transition. The results then depend on the particular location of this transition. 
This drawback can be slightly circumvented by the use of a smooth transition (like ficti- 
tious mixture models in multiphase fiows for instance). Nonetheless, the need to designate 
a specific region where the shift between the two models takes place is detrimental to the 
robustness and the reliability of the model. 

Once the location of the interface or transition region between P'^ and P° has been 
determined, the second problem to solve is the definition of the coupling strategy between 
the two models. Usually, P^ involves some kind of reduction of information compared to 
P^. For instance, in the low-Mach number limit, the velocity field becomes divergence 
free, meaning that it depends on two independent scalar quantities instead of three like in 
compressible fiows. At the interface, in the passage P^ — ?■ P^, it is necessary to project the 
unknowns of P^ onto those of P°, and vice versa, in the passage P° — t- P^ the unknowns of 
P^ must be reconstructed from those of P°. For instance, in the passage from compressible 
to low Mach-number regimes, the velocity must be projected onto a divergence free field. 
In the reverse transition, the irrotational part of the velocity must be reconstructed from 
a divergence free velocity. 

The answer to such questions is by no means obvious and the simulation results also 
depend on the choices of the projection-reconstruction operators. Connection conditions 
can be sought by solving an interior layer problem connecting the state variables of the P^ 
problem at x = — oo to those of the P" problem at x = +00 through a spatial rescaling of 
P^ in the direction normal to the interface. However, quite often, this analysis does not 
lead to a closed set of connection conditions. The interior layer problem itself carries some 
approximations because it involves a rescaling which leads to neglecting all derivatives in 
the tangential direction to the interface. 

Finally, supposing that the questions above have found a satisfactory answer, the whole 
strategy must be practically implemented. The mesh must be constrained to match the 
interface. Therefore, if the interface is not planar, which is very likely in a realistic case, 
an unstructured mesh must be used. Additionally, the interface location may have to 
evolve in time. This brings several additional questions. First, the motion strategy for 
the interface must be defined. What criterion will decide for this motion ? Second, if the 
interface is moved, the mesh must be moved accordingly, otherwise the matching with 
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the interface will be lost. Moving meshes with time is complex, uneasy and costly, both 
in terms of CPU time and of development time. 

All these questions have led an increasing number of teams to look for schemes which 
are valid in both the e = 0(1) and e ^ 1 regimes. For instance, there is an increasing 
literature about so-called 'all-speed schemes' for compressible flows, which are valid for 
all values of the Mach number [71 [25], [30l |45]. However, some precautions must be made 
because, a scheme may be stable in the e <C 1 limit and yet provide a wrong solution, i.e. 
a solution which is not consistent with the P° problem. The correct concept for doing 
so is the so-called Asymptotic-Preserving (AP) scheme, which is described in the next 
section. 



Asymptotic-Preserving methodology 

The concept of an Asymptotic-Preserving method has been introduced by S. Jin [17] for 
transport in diffusive regimes. A scheme P"^'^ for with discretization parameters h 
(standing for both the time and space steps) is called Asymptotic Preserving (or AP) 
if it is stable irrespective of how small the perturbation parameter e is, and if it leads 
to a scheme P^'^ which is consistent with the limit problem P° when e tends to zero 
with fixed discretization parameters h. This property is illustrated in the commutative 
diagram below. 



P' 



e,h fe^O 



> P' 



e->-0 



pO,h 



AP schemes are extremely powerful tools as they allow the use of the same scheme 
to discretize P^ and P°, with fixed discretization parameters. This not possible if the 
AP property is not satisfied because, either the scheme develops an instability when 
the microscopic scale is under-resolved (i.e. when e is too small), or it is stable but 
not consistent with the limit problem P^ and therefore, provides a wrong solution. By 
contrast, when the order of magnitude of e changes dramatically, the AP scheme realizes 
an automatic transition between P^ and P°. There is no model change, no need to 
define a coupling strategy, an interface location, or divising complicated adaptive meshing 
strategies to follow the interface. Thanks to these properties, AP-schemes are extremely 
efficient, producing computational saving of several orders of magnitude. 

The literature about AP schemes is recent, yet increasingly abundant and applied to 
various contexts (see e.g. [SlIllElElIinilMlllSlinilSnilMllSn]). 

In this document, we will illustrate a general methodology for devising AP schemes. 
This methodology is divided into two main steps. The first step is a 'Reformulation' 
step. Indeed, the passage P'^ — >■ P° leads to a change in the type, nature or simply 
expression of the equations which determine some of the unknowns. The reformulation 
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step consists in finding an equivalent form of problem such that explicitly 
appears as a perturbation of P". This concept may sound vague but we hope that the 
various illustrations below will convince the reader of its relevance. To some extent, the 
reformulation procedure consists in bringing the limit model into the model P^ 'by 
brute force'. 

Once the Reformulation step is passed, the Discretization step per se can start. The 
most obvious idea, which is to discretize the reformulated model P*^ should actually be 
discarded. The reason is that P^ is usually quite complicated and it is not clear what 
a 'good' discretization of it is. Additionally, there are some structural properties which 
should be preserved by the discretization. For systems of conservation laws for instance, 
these properties could be a special conservative form or a special formula for the numerical 
viscosities. It is not clear at all how to enforce these properties on a direct discretization 
of R'. 

So, the discretization step proceeds as follows: Discretize (and not P^) into a 
scheme P^''^ in such a way that the various manipulations which led to P"^ in the con- 
tinuous case can be performed at the discrete level. Then, the reformulation procedure, 
performed at the discrete level, permits the derivation of a reformulated scheme P^''* which 
is consistent with the reformulated continuous model P^. Then, the scheme P^''^ appears 
as a perturbation of a scheme R^'^ which is consistent with the limit problem P°, in just 
the same way as the reformulated model P"^ appears as a perturbation of P°. By this 
procedure, in the case of conservation laws, additional properties such as conservativity, 
or special choices of numerical viscosities can be imposed on P'^''* and can be carried over 
to P^''^ by the discrete reformulation step. The discrete reformulation step also leads to 
a form which is the most suited for numerical resolution. 

In other words, it is preferable to 'reformulate the discretization', rather than 'dis- 
cretize the reformulation'. We will illustrate this motto on different examples throughout 
this document. 

We now make some comments. The reformulation procedure requires that the limit 
problem P*^ is well identified and well-posed. In the current state-of-the-art, there is no 
possibility of deriving AP schemes for models P^ whose limit P° is not well-identified and 
well-posed. 

The second remark is that, for the scheme P^'^ to be 'reformulable' into R'^'^ and 
to have a limit P°'^ when e — )■ 0, some level of implicitness is required. The design of 
an AP scheme is just about determining which terms have to be evaluated implicitly in 
order to meet the requirements of the AP property. Sometimes, there is not a unique 
way to perform this goal (we will see an example below). In all cases, this key point is 
the most important and difficult part of the procedure. For this reason, in most of the 
cases, constructing an AP scheme is above-all a question of time discretization. This is 
why, in these notes, we will systematically devise the AP methodology first in the time 
semi-discrete setting, before applying it to the fully-discrete case. 

Then, the question is: why not choosing all terms fully implicit. Then, we are guar- 
anteed to make an AP scheme. There are several issues. First, this is doing too much, 
because requesting a scheme to be AP is not asking for inconditional stability, but only for 
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uniform stability with respect to a given parameter, which is weaker. The fully implicit 
scheme must be inverted. This requires an iterative procedure (because the problem is 
usually nonlinear) inside which a large linear system solver must be called. For systems 
of conservation laws, the linear operator to invert is issued from a first order differen- 
tial system. Therefore, it is generally not positive definite and its inversion can be quite 
unstable. 

By contrast, in designing AP schemes it is usually possible to reduce the inversion 
of the implicit part to the inversion of an elliptic operator. This leads to a much more 
stable numerical resolution. In the case of systems of conservation laws in which the limit 
e ^ corresponds to some characteristic speeds going to infinity, this amounts to solving 
a wave equation (for the fields associated to these characteristic speeds) by a fully implicit 
scheme. Therefore, the computation of the next time iterate involves the inversion of an 
elliptic equation corresponding to the stationary wave equation. 

When dealing with a fully implicit scheme, the iterative procedure needed to invert the 
nonlinear system is usually a Newton method, which is based on a linearization procedure. 
When the limit e — )■ corresponds to the enforcement of some nonlinear constraints, the 
linearization procedure may destroy the AP character of the scheme. This is the case for 
instance in the low Mach-number limit which is characterized by the condition that the 
gradient of the pressure should vanish. For real gases, the pressure is a nonlinear function 
of the state variables. Linearizing this condition by the Newton method may lead to large 
errors. Of course, if full convergence of the Newton iteration is reached, the scheme is 
AP. But this is never the case. The stopping criteria consist in checking if the residual 
of the iterations is below some threshold. But in the limit e — 0, the conditionning of 
the problem is very bad and checking on the residual does not guarantee a reasonable 
convergence level for the solution. Then, the solution is not only bad, but also inconsistent 
with the nonlinear constraint characterizing the limit regime, because of the linearization 
in the Newton iterations. 

We hope that this discussion has been convincing that fully implicit schemes are not 
'the' ultimate method for perturbation problems. Of course, fully implicit scheme may 
be valuable, but they have to be implemented with iterative methods which guarantee 
the AP character of the scheme even if the iterations are stopped before convergence. In 
our experience, going from an AP scheme to a fully implicit scheme only brings minor 
improvements in general. 

Plan of these notes 

In these notes, we plan to develop the AP methodologies on two different classes of 
asymptotic limits. 

Part 1 will be devoted to the quasi-neutral limit in plasma fluid models. We will suc- 
cessively consider the one-fluid isentropic Euler-Poisson model and the one- fluid isentropic 
Euler-Maxwell models. The extension to two or more species, to full Euler models instead 
of isentropic models, and finally, to other variants of the quasi-neutral limit such as the 
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case of the Euler-Poisson-Boltzmann model are briefly discussed at the end. The material 
of these two sections has been (or will be) published in [131 CH EH ESI EEl El] . The appli- 
cation of this strategy to kinetic models (Vlasov-Poisson equations) has been published 
in [21 [231 El] ■ The study of the transition from quasineutrality to non-quasineutrality in 
plasmas has been investigated in [HI [151 EH EH EI] , as well as in [371 [131 [571 [6l] . The 
quasi-neutral limit has been theoretically investigated in [HI [39l [48l [Ml [62] . 

Part 2 will be devoted to the drift-fluid limit in the Euler-Lorentz system, i.e. the 
isentropic Euler system of gas dynamics equations subjected to a Lorentz force. The 
scaling is such that the whole Lorentz force and the pressure force are large compared 
to the inertial force. In this case, we will show that two AP strategies can be divised. 
Both lead to the resolution of a strongly anisotropic diffusion operator corresponding to 
diffusion along the magnetic field lines. The material of this section has been published in 
[20] and is the subject of works in preparation [HI [IT]. Related material for more general 
anisotropic diffusion equations or for the full Euler-Lorentz case is or will be available in 
[I9l [T6l [9l [H]. The drift-fluid limit is the cornerstone of many physical models, see e.g. 
[321 [35l [38l [H] . We refer the reader to the monograph [16] for the physics of magnetic 
plasma confinement. The mathematical investigation of this limit is to the best of our 
knowledge, open. 

In these notes, we focus on the concepts and ideas. We refer the reader to the bibliogra- 
phy given above for applications to practical cases and performance tests. All performance 
tests indicate that these schemes are more efficient than standard schemes by several or- 
ders of magnitude (for instance, a factor 10"^ has been reached in the case of the drift-fluid 
limit, see [20]). 
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Part 1 
Quasineutrality 
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1 The Euler-Poisson system 



1.1 Setting of the problem 

The One-Fluid Euler-Poisson (EP) model describes the plasma electrons through a system 
of isothermal or isentropic gas dynamics equations subjected to an electrostatic force. The 
electrostatic potential is related to the electron density through the Poisson equation. The 
ions are supposed to form a neutralizing background, i.e. they are steady and with uniform 
density. The EP model is written 



Here, n{x,t) > 0, u{x,t) G M.'^ and (f>{x,t) G M stand for the electron density, electron 
velocity and electric potential respectively, and depend on the space-variable x G M'' and 
on the time t > 0. Strictly speaking, the ion density ni{x,t) itself satisfies a system 
of Euler equations. However, in this section, we suppose it uniform and constant in 
time for simplicity and ignore the question of the coupling of the electrons to the ions 
(see section [3]). The positive elementary charge is denoted by e and the electron mass, 
by m. The electron pressure p is supposed to be a given function of the density (e.g. 
p{n) = ksTn in the isothermal case, where T is the electron temperature and ksi the 
Boltzmann constant, or p{n) = CrC , where 7 > 1 and C > are given constants, in the 
isentropic case), eo refers to the vacuum permittivity. The operators V, V- and A are 
respectively the gradient, divergence and Laplace operators and u®u denotes the tensor 
product of the vector u with itself. Finally, d is the dimension {d = 1,2, or 3). 

If the quasi-neutral assumption is made, the Poisson equation (11.31) is replaced by the 
constraint of zero local charge: 



In this context, the quasi-neutral Euler-Poisson model coincides with the incompressible 
Euler (IE) model: 



together with (II. 4p . with the hydrostatic pressure tt = —em~^(j). In deriving (II. 6p . we 
have used that for u satisfying (11.51) . V(u ® u) = (m ■ V)n. 

The passage from EP to IE can be understood by a suitable scaling of the model, 
which highlights the role of the scaled Debye length: 



dtu + V ■ (nu) 

m{dt{nu) + V 
g 

-A0 = — (ni ■ 




(1.1) 
(1.2) 

(1.3) 



n = rii 



(1.4) 



V-n = 0, 

dtu + (n ■ V)n + Vtt = 0, 



(1.5) 
(1.6) 




(1.7) 
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where xq is the typical size of the system under consideration. Xd measures the spatial 
scale associated with the electrostatic interaction between the particles. The dimensionless 
parameter A is usually small, which formalizes the fact that the electrostatic interaction 
occurs at spatial scales which are much smaller than the usual scales of interest. However, 
there are situations, for instance in boundary layers, or at the plasma-vacuum interface, 
where the electrostatic interaction scale must be taken into account. This means that 
the choice of the relevant macroscopic length xq may depend on the location inside the 
system and that in general, the parameter A may vary by orders of magnitude from one 
part of the domain to another one. 



1.2 Scaling of the EP model and quasineutral limit 

Let Xq, to, Wo, Po, 00 and uq be space, time, velocity, pressure, potential and density scales. 
Scaled position, time, velocity, pressure, potential and density are defined hj x = x/xq, 
t = t/to, u = u/uq, p{n) = p{n)/pQ, cj) = —4>/4>o and n = n/uQ. We choose Xq to be 
the typical size of the system (for instance an inter-electrode distance or the size of the 
vacuum chamber). We also choose riQ = n^. We define a temperature scale by the relation 
Pq = uokBTo and define the velocity scale as uq = {kBTom"^)^/'^ . Finally, the potential 
scale is set to 0o = ksTo/e, the so-called thermal potential. 

Inserting this scaling and omitting the bars gives rise to the EP model in scaled form: 

Definition 1.1 The scaled Euler-Poisson (EP) model is: 

+ V ■ (nV) = 0, (1.8) 
dtin^u^) + V ■ {n\^ ® u^) + Vp{n^) = n^V4>^, (1.9) 
-\^A(f)^ = l~n^. (1.10) 



where A is the scaled Debye length ( 1.1 ). 



Formally passing to the limit A — )■ in this system and supposing that — )■ n*^, m'^ — 
(f)^ — >■ 0° as smoothly as needed, we are led to the Icompressible Euler (IE) model: 

Definition 1.2 The scaled Incompressible Euler (IE) model is: 

V-M° = 0, (1.11) 
+ (m° ■ V)m° = V0°, (1.12) 
n° = 1. (1.13) 



Indeed, the first two equations are the scaled form of the Icompressible Euler (IE) model 
where the potential cjP is the negative hydrostatic pressure and acts as the Lagrange 
multiplier of the incompressibility constraint (11.111) . In order to resolve the constraint 
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and find an equation for 0", we take tlie divergence of (11.121) and insert tlie constraint 
(fLTTD . We find: 

A0° = V' : (m°®m°). (1.14) 

Tliis equation is equivalent to (11.111) provided tliat (II. lip is satisfied initially. Therefore, 
we state: 

Definition 1.3 The Reformulated Incompressible Euler (RIE) model is: 

n° = 1. (1.15) 

+ (m° ■ V)m° = V0°, (1.16) 

A0° = V': (m°®m°). (1.17) 



Proposition 1.4 The RIE model supplemented with eq. 1^1. at time t = is equivalent 
to the IE model. 

The question arises if tliere exists a Reformulated Euler-Poisson model in the same 
way as there exists a Reformulated Incompressible Euler model. Using the mass and 
momentum conservation equations, such a reformulation can be derived. Indeed, taking 
the time derivative of (11. Sp . the divergence of ( 11. 9p and eliminating the second order time 
derivative of n'^ by means of the Poisson equation (ll.lOp . we are led to the so-called 
'Reformulated Poisson equation': 

-V ■ [{n^ + X'^d'^)V(j)^] = -V^ : {n\^ ® + p{n^)ld) , (1.18) 

We refer to the Reformulated Euler-Poisson (REP) model as the model consisting of the 
mass and momentum conservation equations complemented with the reformulated Poisson 
equation (11.180 : 

Definition 1.5 The Reformulated Euler-Poisson (REP) model is: 

^tn^ + V ■ (n V) = 0, (1.19) 
dt{n\^) + V ■ {n\^ ® u^) + d^p{n^) = n^V<p^, (1.20) 
-V ■ [{n^ + X^d^t)V<i)^] = -V^ : [n^u^ ® + p{n^)Id) . (1.21) 

A solution of the REP model is a solution of the original EP model if and only if the 
initial data satisfy the Poisson equation in its original form (ll.lOp together with its time 
derivative. Indeed, (11.18P is a second order differential equation in time and requires 
Cauchy data for 0° and its time-derivative at t = 0. The equivalence between the EP and 
REP models is summarized in: 
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Proposition 1.6 The REP model is equivalent to the original EP model ( f j.<§l) . ( f i.gj) . 
( li. iOj) . provided that the Poisson equation in its original form l\1.10\) and its time deriva- 
tive are satisfied at initial time. 

Note however that, at the discrete level, schemes based on the REP model may differ from 
those based on the EP model. It clearly appears, that f ll.2ip is a singular perturbation 
of (11 -Mp when A — )■ and consequently that the REP model formally tends to the RIE 
model in this limit. Therefore, using the REP formulation appears as a good strategy 
to devise schemes for the EP model which are consistent with the IE model in the limit 
A — )■ 0. In the next section, we propose a time semi-discrete scheme based on this strategy. 

1.3 Time-semi-discretization and AP property 

We denote by 6 the time step. For any function g{x,t), we denote by g"''{x) an approx- 
imation of g{x,t"^) at time t"^ = m6. The classical time-semi-discretization is based on 
the EP model. We propose a new time-semi-discretization based on the REP model. 

In the classical method, the force term in the momentum equation is taken implicitly. 
S. Fabre [M] has shown that this implicitness is needed for the stability of the scheme 
(an explicit treatment of the force term leads to an inconditionally unstable scheme). 
Additionally, this implicitness still gives rise to an explicit resolution, since the mass 
conservation can be used to update the density, then the Poisson equation is used to 
update the potential, and finally the resulting potential is inserted in the momentum 
equation to update the velocity. However, stability is subjected to a condition of the type 
S = 0(A) and the method breaks down in the quasineutral limit A — )■ 0. 

Definition 1.7 The classical time- semi- discretization of the EP model is: 

^-l(^A,m+l _ ^A,m^ ^ y . (n^'"^M^''") = 0, (1.22) 
^-l(^A,m+l^A,m+l _ ^A,m^A,m^ ^ y . (^^''^M^''" ® M^'™) + Vp(n^'"^) = 

= (1.23) 

-A^A^^'^+i = 1 - n^'^+i. (1.24) 



This scheme cannot be AP. Indeed, the limit A — )■ with fixed 5 leads to 

^_1^^0,m+l _ ^0,m^) ^ y . (riO'-u"''") = 0, (1.25) 

= n°''^+V0°''"+\ (1.26) 
= 1, (1.27) 

which does not provide a valid recursion for computing the unknowns at time t"^~^^ knowing 
them at time t"^. Indeed, (ll.27p . once inserted into fll.25p . provides a divergence constraint 
on M^'"^, which is impossible to fulfill since is a datum from the previous time step. 

A cure for this deficiency is to evaluate the mass flux implicitly. Simultaneously, we 
can see that the density in the force term can be taken explicit. This leads to the following 
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Definition 1.8 The AP time Semi-Discretization of the EP model (or SD-EP scheme) 
is: 

^-l^^A,m+l _ ^\m^ ^ y . (^A,m+l^A,m+l _ ^^ 28) 

= n^''^V0^''"+\ (1.29) 

-A2A(/.^'™+i = 1 - (1.30) 

Using the discrete momentum eq. f ll.29p and the Poisson eq. (11.301) respectively, n^^^^^ 
^A,m+i ^A,m+i ^^^^ eliminated from the mass equation (11.281) . This leads to: 

-V ■ ((^^n^''" + A2)V(/)^''"+^) = 1 - n^'" + 5V ■ (n^''"u^''") 

-5^V ■ (V ■ (n^'™M^'™ ® M^'™) + Vp(n^''")) . (1.31) 

This elliptic equation for is a discrete analog of the reformulated Poisson eq. (I1.2ip . 

Indeed, using (I1.28P between times m — 1 and m and (11.301) . we can write it as 

-V ■ ((5^^^'™ + A2)V0^'™+^) = A2(-2A0^'™ + A0^''^-i) 

-^^V^ : ® u^''" + p{n^'"')ld) , (1.32) 

where the discretization of the second order time-derivative of A0 appears. We summarize 
this in the 

Definition 1.9 The Reformulated time Semi- Discretization of the EP model (RSD-EP 
scheme) is: 

^_l^^A,m+l _ ^A,m^ ^ y . ^^A,m+l^A,m+l ^ _ ^^ 33) 

^_l^^A,m+l^A,m+l _ ^A.m^A.m^j ^ y . (^^A,m^A,m ^ ^A,m^) ^ Vp(n^''") = 

= n^'"^V0^'™+\ (1.34) 

-V ■ ((<52n^'™ + A2)V(/>^'"'+^) = 1 - n^'™ + 5V ■ (n^'"^u^'"^) 

-5^ ■ (V ■ (n^'^M^'"* ® M^'™) + Vp(n^'™)) . (1.35) 



Proposition 1.10 The RSD-EP scheme is equivalent to the SD-EP scheme. 

Formally passing to the limit A — )■ with fixed 5 in this scheme leads to the following 
scheme: 
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Proposition 1.11 We consider a sequence of solutions of the RSD-EP scheme whose 
initial conditions are well-prepared, i.e. satisfy n^'^ — 1 = 0{\^), V ■ {n^'^u^'^) = 0{\^). 
Then, the formal X ^ of this sequence of solutions satisfies the following scheme: 

= 1, (1.36) 

^-l(^0,m+l _ ^0,rn^ ^ y . (^0,m ^ ^0,m^ _ V0°'™+\ (1.37) 

r/iis scheme (the Reformulated Semi-discrete Incompressible Euler scheme or RSD-IE 
scheme) is consistent with the IE model. 

Proof: Leting A — )■ in the RSD-EP scheme leads to: 

^_1^^0,m+l _ ^0,m) ^ y . ^^0,m+1^0,m+l) _ g, (1.39) 

= n°'"^V0°''"+\ (1.40) 

-5^ ■ (n°'™V0°''"+^) = 1 - ri°'™ + 5V ■ (n°''^M°'"*) 

-<52V • (V • (n^'^n^-™ ® u°'™) + Vp(n°''")) . (1.41) 

Taking the combination 5 x ffL39|) - 5^ x V ■ fOOD + fOTD leads to f06|) for m > 0. 
With the well-prepared initial condition assumption, we have n^'^ = 1 for all m > 
and consequently, using (11.391) again, to V ■ (t^O'^+i^^O'^+I) = Q for all m > 0. With the 
well-prepared initial condition again, we deduce that V ■ (T^O'^-yO.'^) = Q for all m > 0. 
Then, (jLiTjl reduces to (|L38|1 for m > 0, while (iLinil reduces to (fOTIl for m > 0. This 
shows that the scheme is consistent with the RIE model. ■ 



Remark 1.1 The well-prepared initial condition hypothesis on n is satisfied as soon as 
the first iterate m = satisfies the Poisson eq. —X^^cj)^'^ = 1 — n^'^ which is a natural 
condition for a solution of the Euler-Poisson problem. If the well-prepared initial condition 
hypothesis on u^'^ is not satisfied, the RSD-IE scheme is still satisfied but starting at iterate 
m = 1. If none of the well-prepared initial condition hypotheses is satisfied, the RSD-IE 
scheme is satisfied starting at iterate m = 2. 

This proposition and the previous ones show that, in the A — ?■ limit, the SD-EP scheme 
is consistent with the Incompressible Euler model. Therefore, the SD-EP is AP, provided 
that its stability condition is independent of A when A — )■ (Asymptotic Stability prop- 
erty). This will be shown, at least in a linearized setting, by the stability analysis below. 
We note that the computational cost of the AP-scheme is the same as for the classical 
scheme, the Poisson equation being replaced by the elliptic equation (ll.3ip . which has 
roughly the same computational cost. This AP scheme has been proposed initially in [13] 
and [11]. 
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1.4 Linearized stability analysis of the time-semi-discretization 

The linearized EP model will be considered about the state defined hj = 1, = 0, 
(1)^ = (which is obviously a stationary solution). We will also restrict ourselves to a 
one-dimensional situation. Expanding n'*' = 1 + en^, = eu^ and cf)^ = e(j)^, with e 
being the intensity of the perturbation to the stationary state, and retaining only the 
linear terms in e, we find the linearized EP model: 

dth^ + d^u^ = 0, (1.42) 
dtu^ + Td,fi^ = d,^\ (1.43) 
X^dl4>^ = n\ (1.44) 

with T = p'{l). Introducing n^, u^, (p^, the partial Fourier transforms of n^, u^, cf)^ with 
respect to x, we are led to the system of ODE's: 

dth^ + i^u^ = (1.45) 

dtu^ + iTin^ = 1^4)^, (1.46) 
-A^^^^A _ ^A^ ^^ 4-^^ 

where ^ is the Fourier dual variable to x. The general solution of this model takes the 
form 



with 

4 = ±z(Te2 + A-2), (1.49) 

and n^, u\ are given functions of ^, fixed by the initial conditions of the problem. In 
particular, since s\. are pure imaginary numbers, the norm of the solution is preserved 
with time. When A — > 0, then |s^| — )■ oo, which means that the limit A — !■ is an 
oscillatory limit. We will see that the AP-scheme averages out these fast oscillations. 

The goal of this section is to analyze the linearized stability properties of both the 
classical and SD-EP schemes. We want to show that the classical scheme requires a CFL 
condition of the type 5 = 0(A) while the stability condition for the SD-EP scheme is 
independent of A when A — )■ 0. This last property is known as 'Asymptotic-Stability' 
and is a component of the Asymptotic-Preserving property (see introduction). Indeed, 
the faculty of letting A — )■ in the scheme with fixed S is possible only if the stability 
condition of the scheme is independent of A in this limit. We will prove L^-stability 
uniformly with respect to A for the linearized problem fll.45l) - fll.47l) . 

In general, time-semi-discretizations of hyperbolic problems are inconditionally unsta- 
ble. This is easily verified on the discretization fll.45p - fll.46p if is set to 0. This is because 
the skew adjoint operator dx has the same effect as a centered space-differencing. For fully 
discrete schemes, stability is obtained at the price of decentering, i.e. adding numerical 
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viscosity. To mimic the effect of tliis viscosity, in the present section, we will consider 
the linearized Viscous Euler-Poisson (VEP) model, which consists of the linearized EP 
model fll.42p - fll.44p with additional viscosity terms (in this section, we drop the tildes for 
notational convenience): 

dtn^ + d^u^ - /3dln^ = 0, (1.50) 
dtu^ + Td^fi^ - Pdlu^ = d^cf)^, (1.51) 
A'a^0^ = n\ (1.52) 

where /3 is a numerical viscosity coefficient. We keep in mind that, in the spatially 
discretized case, /3 proportional to the mesh size h: 

13 = ch, (1.53) 

with the constant c to be specified later on. 

The classical and SD-EP time discretizations are given by 

^-i(^A,m+i _ ^x,m^ ^ - pdln^'"" = 0, (1.54) 

^-i^^x,m+i _ ^x,m^ ^ Td.n^'"' - (3dlu^''^ = (1.55) 

^2Q2jX,m+l ^^X,m+1^ (1.56) 

with a = for the classical scheme and a = 1 for the SD-EP scheme. Passing to Fourier 
space with C, being the dual variable to x, and eliminating 0'*'''""'"^, we find the following 
recursion relations: 

A,m+1 _ ^A,m^ ^ i^u^'"'+^ + /J^V''" = 0, (1.57) 

^-i^^A,m+i _ ^A,m) ^ iT^n^^"" + ^n^'"+' + /^e'u^'"' = 0, (1.58) 
The characteristic equations for thus recursion formula is: 

-2q(^l- (^eS - ^) + (1 - + TeS' = 0. (1.59) 

for the classical scheme (a = 0) and 



A 



- 2g ( 1 - pes - ;;Tes' ) + (1 - pesy = o. (i.eo) 



for the SD-EP scheme (a = 1), where q is the characteristic root. Each of these quadratic 
equations has two roots q^{C,) which provide the two independent solutions of the corre- 
sponding recursion formulas. Their most general solution is of the form 

( £11) ^ ( l«) ) ■ 
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where n^{C,) and u^{^) depend on the initial condition only. A necessary and sufficient 
condition for stability is that < 1- However, requesting this condition for all 

^ e M is too restrictive. To account for the effect of a spatial discretization in this analysis, 
we must restrict the range of admissible Fourier wave- vectors to the interval [— f , f]- 
Indeed, a space discretization of step h cannot represent wave- vectors of magnitude larger 
than ^. This motivates the following definition of stability: 

Definition 1.12 The scheme is stable if and only if 

\qiiO\<^^ such that 1^1 < ^- (1-62) 

Now, our goal is to find sufficient conditions on S such that either schemes are stable. 
More precisely. In [26], the following result is proven: 

Proposition 1.13 We assume that the numerical viscosity fi satisfies 
(i) We assume that 5 satisfies the CFL condition 5 < CqH with Cq = jr^:^. Then, there 
exist two constants < Ci < C[ and a constant Ai > (depending on c, T), such that 
the classical scheme is stable provided that 6 < CiX and unstable if S > C[\, for all \, 
< A < Ai. 

(a) There exists a constant C2 > (depending on c, and T but not on X), such that the 
SD-EP scheme is stable for all X > provided that 6 < C2h. 

The ffist statement of this theorem says that, under the CFL condition for the hydro- 
dynamic part of the system, a necessary and sufficient condition for the stability of 
the classical scheme is that 6 < 0(A). This conffims that the classical scheme is not 
Asymptotically-Stable and cannot be AP. 

The second statement shows that the stability condition of the SD-EP scheme is 
independent of A in the limit A — (and actually, whatever value A has). Thus the 
Asymptotic Stability property is fulfilled. The stability is not unconditional: a CFL 
condition for the hydrodynamic part is still required. However, the goal of an AP scheme 
is not to provide unconditional stability, but simply Asymptotic Stability with respect to 
a given parameter. This more restricted stablity requirement allows us to select the terms 
that need to be evaluated implicitly, and leads to schemes which are more easily resolved. 

Adding more implicitness in the scheme does not necessarily improve the AP-ness 
of the scheme, and sometimes a degradation is observed. Indeed, fully implicit schemes 
require nonlinear solvers which furnish only an approximate solution of the scheme. They 
may perform poorly in highly ill-conditionned situations; The stopping test in case of 
iterative solvers may also be not as accurate as expected. The result is that the AP 
property can actually be missed by fully implicit solvers. It can be retrieved at the price 
of some specific AP pre-conditionning. 

The linearized stability analysis does not provide a complete proof of the AP property. 
For this purpose, energy estimates on the fully non-linear system should be proved. Some 
partial answers can be found in [2H]. The stability analysis can also be developed about 
a stationary state with non-zero velocity, with a similar conclusion. We refer the reader 
to |26] for more detail. 
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1.5 Full (space-time) discretization: enforcing the Gauss law 

1.5.1 General framework and classical fully-discrete scheme 

We now discuss some issues related to tlie spatial discretisation. For this purpose, we 
restrict ourselves to a one-dimensional framework but all the considerations here extend 
straightforwardly to any dimension. We introduce a spatial discretization with a uniform 
mesh of step h and we denote by Ck the cell [{k — l/2)h, {k + l/2)h] and Xk = kh, with 
/c G Z. Like in usual first-order shock capturing schemes, the fluid unknowns n and u are 
approximated by piecewise constant functions within the cell Ck and represented by cell- 
centered values n|^, at time = mS. The electric potential is also approximated 
by cell-centered values The discretization of the hydrodynamic part is performed by 
means of a flrst order shock capturing scheme. We denote by /n 1^^1/2 /"lfc+1/2 

the 

numerical fluxes for the mass and momentum conservation equations respectively, at time 
_ ^YiQ cell interface Xk+i/2- Similarly, -E|^^i/2 denotes the approximation of 

the electric fleld E = —d^cf) at Xk+1/2 and t"^. 
The classical scheme is deflned as follows: 

Definition 1.14 The fully- discretized classical scheme for the Euler-Poisson model is: 

5~\n\^+^ -n\^) + /^-^(/„|r+i/2 - /niri/2) = 0, (1.63) 
ri((n«)ir' - + h~\fX^y, - fX-1,2) = (1.64) 

X'h-\E\^^l,, - £;|-Y/2) = 1 - (1-65) 

^ir+v2 = -h-\K:i - ^ir' = ^(^ir;i/2 + ^1^1/2)- (i-66) 

The numerical fluxes fn\'k+i/2 ^'^^ /«lfc+i/2 ^'^^ given by: 

/nlr+1/2 ^ 1 r {nu)\^ \ ( (nn)l-, 



+ 



( rnM)|^-(nt)V+i )^ ^^-^^^ 



where /x^;^ ^ is the viscosity matrix. 



Here, we will not discuss the choice of the viscosity matrix (see jlQl [531 ES] fo'^ a discussion 
of this point). For instance, in the simple case of a Local Lax-Friedrichs (LLF) scheme 
[531 158] corresponds to a local evaluation of the wave speeds (see e.g. [H] for more 

details). 

We note that the electric fleld is implicit in the momentum conservation eq. As 
shown by Fabre [M], this is a necessary condition for conditional stability. In spite of 
this implicitness, the scheme can be solved explicitly. The sequence of updates goes as 
follows: flrst the mass conservation eq. f ll.63p is used to compute n™'"'"^. Then (I1.65P and 
(11.661) are combined into the discrete Poisson equation 
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the inversion of which allows us to compute (provided suitable boundary conditions 
are specified. We will not discuss this point here). Finally, with the momentum balance 
eq. f ll.64p we find the value of {nu)\'^~^^. The stability condition for this scheme combines 
a CFL condition for the hydrodynamic part and a condition of the type 6 < 0(A) [31] . 
Therefore, this scheme cannot be AP. 

1.5.2 AP fully discrete scheme 

The AP-scheme is defined as follows: 

Definition 1.15 The Fully Discretized AP scheme for the Euler-Poisson model (FD-EP 
scheme) is: 

5~\n\r' - + - /nirV2) = 0, (1.69) 

5-\{nu)\^^' - {nu)\^^) + /i-^(Mr+i/2 - M^-i/s) = -n\tE\T\ (1-70) 

= -h-\K:r' - <t>\r% ^ir^ = \{e\^^i„ + ^^in^)- (1-72) 

The momentum flux fu\"'k+ii2 '^^ given by jl-dl^ while the mass flux is given by: 



Jn\k+l/2 ~ ■l'n\k+l/2 2 ^ 2 '^•'"l*''+3/2 Ju\k-l/2)- \^-''^J 



Compared to the classical scheme, the FD-EP scheme exhibits a modified mass fiux fll.73p . 
and uses an explicit evaluation of the density appearing in front of the electric field in the 
momentum equation (11.701) . We now provide the rationale for the use of the mass fiux 
(11.731) . In the following discussion, we will restrict ourselves to an LLF (or Rusanov) scalar 
numerical viscosity for simplicity, but the discussion could be extended straightforwardly 
to any kind of shock-capturing method. We leave the details to the reader. We show that 

the fiux (11.731) can be interpreted as some minor modification of a fiux fn\Z+i/2 which 
the central part is implicit: 

/nir+Y/2 = I [inu.)\r' + inu.)\T:i + (^ir - ^ir+o] • (1-74) 

Indeed, starting from the definition (ll.74p and inserting (11.701) into (11.741) . we can write: 

6 



~ Jn\k+l/2 ~ ^ 



I fc+1 -£^1^+3/2 + I'^lfc+l + "-Ifc J -^lfc+1/2 + -^lfc-1/2 

7y (/w|fe+3/2 ~ Ju\k-\I2)- (1-75) 



This fiux involves an average of £""+1 over three neighbouring mesh points. In order to 
reduce the numerical diffusion, we replace (11.750 by (ll.73p . Therefore, (11.730 is, up to 
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a lumping of the electric field average, what results from an implicit evaluation of the 
central part of the mass flux. It is an order 0{6) modification of the explicit flux, but 
this simple modification is crucial in making the scheme AP. Formula (I1.73P can be used 
irrespective of the expression of the explicit fluxes /n|^i/2' /«lfc+i/2 therefore valid 

whatever choice of numerical viscosity is made. 

We now show that the FD-EP scheme can be reformulated in a consistent way to the 
REP model. 



1.5.3 Reformulation 

Definition 1.16 The Reformulated Fully Discretized AP scheme for the Euler-Poisson 
model (RFD-EP scheme) is: 



fu 



0, 

m \ 
k-1/2) 



-n\TE\r\ 



-jin\T+i + n\m<P\ 



m+l 
fc+1 



|m+l^ 



52 

2 



-(A' + ^('^ir+'^ri))('/'ir' 



|m+l^ 



n\ 



h-'5{f„ 



fc+1/2 



fn 



\k-l/2) 



piim+l 
-^lfc+1/2 



im+l^ 
\k 



■U«lfc+3/2 Ju\k+l/2 Ju\k-l/2 "T Ju\k-3/2)- 



^|m+l 



2 V-^lfc+1/2 ^ ^\k-l/2)^ 



;i.76) 
;i.77) 



;i.78) 
;i.79) 



with /n|fe!f4/2 ^''T'd fu\^^i/2 Qivcu by ^1.73\ ) and by ( [i.g7| ) respectively. 



Proposition 1.17 The FD-EP and RFD-EP schemes are equivalent. 
Proof: We first insert (11.691) into (11.711) and get 



im+l 



(A E\k+l/2 '^/"lfc+1/2^ 

Using (ll.73p . we have: 

^2 



(A ^ir-1/2 - 



im+l \ 
lfc-1/2/ 



\2r-i|m+l I m+l 

-^lfc+1/2 "jn\k+l/2 



(A^ 



(n 



-6 



fc+1 ''■Ifc ))^\k+l/2 

6h-^ 



fn 



+3/2 

Then, inserting ( ILSTD and (fL72D into (0(1 leads to ( ILTSD . 



Ju\k-1 



12) 



;i.8o) 



:i.8r 
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Eq. f ll.78p is a discrete elliptic equation which allows to find 0"*"'"^, provided boundary 
conditions (the same as for the Poisson equation) are specified. It clearly appears as a 
consistent spatial discretization of f ll.3ip . The time update for the AP-scheme follows a 
different sequence from that of the classical scheme. We first solve for 0™+^ by inverting 
the discrete elliptic equation fll.78p . Then, u"^~^^ and n™+^ can be explicitly computed by 
successively using (I1.77P and (]1.76p . 

Although not explicitly solved, the Gauss equation is satisfied exactly. More precisely, 
we have: 



Proposition 1.18 The solution of the RFD-EP scheme satisfies the discrete Gauss eq. 
^l.'71\) exactly. 

Proof: It follows the steps of the proof of proposition II . 1 71 backwards. Inserting the first 
equation f ll.79p into f ll.78p . and having f ll.73p in mind, we see that the latter is equivalent 
to 

Inserting the discrete mass balance eq. (I1.76P into f ll.82p leads to the discrete Gauss 
equation f ll.7ip . ■ 

We now investigate the limit A -> 0. We define: 



1.5.4 A —7- limit and AP property 

Proposition 1.19 (i) We consider a sequence of solutions of the RFD-EP scheme whose 
initial conditions are well-prepared, i.e. satisfy n^'^ — 1 = O(A^). Then, the formal X — )■ 
of this sequence of solutions satisfies the following scheme: 



„ im+l 
"-Ik 



m+1 
k 



U 



+ h-\fu\ 
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I m+1 
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(1.84) 
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;i.87) 



/^fc+i/2 ^^mg' the viscosity matrix. 
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(ii) We denote by: 



D \ / n 



J-^ll 7^ I 1 /O / 



Then \1.83\) is equivalent to 

-2ml m 1 m 1 m 

— h ^ — 201™"*""^ + 01™.^!"^) = ^(/«lfc+3/2 ~ /w|fc+l/2 ~ fu\'k-l/2 + /^l™- 

(01- , - 201- + 01™ ,) + ^ (01-2 - 20ir + 01^2) 



■3/2J 



-At l-^nU.+l/2 -^"lfe-l/2 -^^"1^+1/2 T" -'^n|fe-l/2J' l-"-' 



for m > 1. // i/ie additional 'well-prepared' initial condition hypothesis h ^1 



u 



k+l 



U 



= O(A^) made, eq. U.88\) holds true for m = without the last two terms. 
Then, this scheme (the Reformulated Fully- Discrete Incompressible Euler scheme or RFD- 
lE scheme ) is consistent with the RIE model, as soon as there exist two constants Ci and 
C2 such that and Cih < S < C2h. If the viscosity matrix is such that E)n\^^i/2 ~ ^ /^'^ 
k, m, then the RFD-IE scheme is consistent wit the RIE model without any condition on 
the time step. 

Proof: (i) Taking the limit A — )■ in the RFD-EP scheme leads to: 

5~\n\^^^ - n|n + /^~'(/nir+Y/2 - /nirv2) = 0' (1-89) 
5~\{nu)\t^' - MID + /^"^(/Jr+i/2 - /Jr-1/2) = -^ir^ir', (1-90) 

= 1 ~ '^Yk + ^ ^'^(/nlfc+1/2 ^ fn\k-~l/2) 

fc+3/2 fu\k+l/2 fu\k~l/2 + fu\k^Z/2)- (1-91) 

^ir+v2 = -^"^(011^1^ - <p\T\ ^ir' = \{E\t:i/2 + ^iri/2)> (1-92) 

We now perform the same calculations as for the proof of proposition 11.181 Inserting 
the first equation (11.921) into (11.911) . and having (I1.73P in mind, we see that the latter is 
equivalent to 

= 1 - + h-H{h\^^^^, - (1.93) 

Inserting the discrete mass balance eq. (I1.89P into (11.931) leads to (ll.83p . With the well- 
prepared initial condition assumption, we deduce that = \ for all m > 0. Inserting 
this information into (ll.90p and (11.670 leads to (11.840 and (ll.87p . For the momentum 
flux, we have removed the constant p(l) because only differences of fluxes at successive 
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mesh points are used and the constant p{l) is cancelled in this process. Then, for m > 1, 
(11.911) implies fll.SSp . This ends the proof of point (i). 

(ii) Taking half the difference of f ll.84p for k + 1 and /c — 1, we obtain: 

r— 1// /• im+l r im+l \ / r \m r \m \\ 

U/n|fe+l/2 \Jn\k+l/2 Jn\k-l/2)) 

r 1 

( p) im+l pi im+l N ^'n I™ _ n \ I 

2 [y^n\k_^.i/2 ^n\k-l/2l \^n\k+l/2 -^nlfc-l/2^J 

"T 2 VWmIA:+3/2 Ju\k+l/2) \Ju\k-l/2 Ju\k-3/2)) ~ 

= ^ (01^+^2' - 20ir^ + 01^2^) , (1.94) 
Multiplying this equation by h^^ and subtracting fll.SSp to it leads to 

s-'h-\fn\T:i)2 - /niri/2) = -h-' {<i>\t:i - 20ir^ + ^/-irY) + 

- 20ir ' + 01^2') 

|m+l p) im+l N ^'n I™' n l^^+l "l 

lfc+1/2 -^"lfc-1/2/ l-^nlfc+1/2 '^n\k^l/2) 



f im+l \ 
in|fc-l/2/ 












+ 2 



;i.95) 



for all m > 0. Inserting fll.QSp into f ll.SSp leads to f ll.SSp for m > 1. For m = 0, the well 
prepared initial condition hypothesis on the velocity leads to the same eq. f ll.SSp without 
the last three terms. Now, the second and third terms of f ll.SSp are 0{h^) approximations 
of dl(p{xk,t"^)- Therefore, their difference is 0{h^). 

We now turn to the last term at the right-hand side of fll.SSp . If the viscosity matrix is 
such that -D„|^^^2 = for all k, m, then, this term is identically zero and what precedes 
shows that fll.SSp is a consistent approximation of fll.l7p . If this is not the case, we must 
estimate the last term at the right-hand side of fll.SSp . For this purpose, let us denote by 
Dn{x,t) a function interpolating -D„|^-^^2 time and space. Then, Taylor's expansion 
up to second order about (xfc,t™~^^/^) shows that 
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Oih{h' + 5')), (1.96) 



because the numerical viscosity is 0{h). Then, the third term at the right-hand side of 
( I1.88P is 0{6 + h'^6~^), which is 0{h) if there exist two constants Ci and C2 such that and 
Cih < 6 < C2h. Under this condition, fll.SSp is still a consistent approximation of fll.lTp . 

That eqs. fll.831) . fll.84p are consistent approximations of fll.lSp and (11.161) respectively 
is obvious. We conclude that the RFD-IE scheme is consistent with the IE model. This 
ends the proof of the proposition. ■ 

The previous proposition shows that the Fully-Discrete FD-EP scheme is AP. The as- 
sumption that Dn\^_^i/2 ~ ^' '"^j satisfied if the viscosity matrix is a scalar, like 
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in the Rusanov scheme, or more generally, if the artificial viscosity applying to the density 
equation only involves the density. For instance, taking any arbitrary viscosity matrix, 
and replacing the line corresponding to the density by only a diagonal term would also 
satisfy this requirement. If this condition is not satisfied, then the AP property requires 
some "inverse CFL condition" Cih < 6 to hold. 

1.6 Euler-Poisson problem: conclusion 

In this section, we have provided an Asymptotic-Preserving (AP) scheme for the one-fluid 
Euler-Poisson problem in the quasineutral limit, i.e. when the scaled Debye length tends 
to zero. In this limit, the Euler-Poisson problem reduces to the Incompressible Euler 
problem. To construct such a scheme, the basic idea is to derive a reformulation of the 
Poisson equation which is not singular when A — > 0. The AP scheme is therefore designed 
to mimic this reformulated Poisson equation. This is done by conveniently choosing the 
terms that must be evaluated implicitly. The AP scheme can be based on any classical 
shock capturing scheme. It consists in perturbing the standard hydrodynamic fluxes by 
corrective terms which are of the order of the time step or the mesh step, and which 
therefore do not emperish the consistency of the scheme. But in the limit A — > 0, they 
do provide consistent approximations of the Incompressible Euler equations. A linearized 
stability analysis has been reviewed. It confirms that the stability condition of the AP- 
scheme is independent of A when A — ?■ 0. However, a nonlinear analysis of the stability of 
the scheme is still lacking. 

2 The Euler-Maxwell system 

2.1 The Euler-Maxwell system and its quasi-neutral limit 
2.1.1 The scaled Euler-Maxwell system 

The one-fluid Euler-Maxwell (EM) system consists of the mass and momentum balance 
equations for the electron fluid coupled to the Maxwell equations. For simplicity, we state 
the problem in an already scaled form. Let n^{x,t) > and u^{x, t) e M"* stand for the 
electron density and electron velocity respectively, where A > is the scaled Debye length 
(11. 7p . They depend on the space- variable x G M"^ and on the time t > 0. The electron 
pressure p = p{n) is supposed to be a given function of n (isentropic assumption) for 
simplicity. We assume that the dimension d = 3 for this presentation. The electric field 
E^{x,t) G R'^ and the magnetic field B'^{x, t) G M.'^ are the two components of the electro- 
magnetic field. The electric sources are the electrical charge p^{x,t) G M and the electric 
current j'^{x,t) G W^. We assume a neutralizing background of immobile ions of constant 
density equal to 1 in scaled units, and of average velocity 0. With these definitions, the 
scaled Euler-Maxwell model is defined as follows: 
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Definition 2.1 The scaled Euler- Maxwell (EM) model is: 

dtn^ + V ■ (n^M^) = 0, (2.1) 

9t(nV) + V ■ (nV ® u^) + Vp(n^) = -n^{E^ + x B^), (2.2) 

dtB^ + V X = 0, (2.3) 

X'^dtE^ -V X B^=f := n V, (2.4) 

V-5^ = 0, (2.5) 

• = := 1 - n^, (2.6) 



Eqs (12. ip . (12. 2p are the mass and momentum balance equations for the electrons. The 
right-hand side of (12. 2 p is the Lorentz force, which depends on the electro-magnetic field 
{E^,B^). The latter is a solution of the Maxwell equations (Q-dM]). Classically, f l23|) . 
(I2.4p are respecctively the Faraday and Ampere equations, and provide the time evolution 
of the electro-magnetic field. The Gauss and Ampere equations (12. 6p and (12. 4p respec- 
tively involve the electric sources {p^,j^) which depend on the hydrodynamical quantities 

and u'^. A lot more physics could be considered. For instance, instead of an isen- 
tropic gas equation-of-state, we could consider a full hydrodynamic model including an 
evolution equation for the gas total energy. We have also neglected electron-ion collisions 
which otherwise would introduce a friction term in (12.21) . Interactions with a neutral gas 
component could also be introduced. The present setting retains the features which will 
be important for the forthcoming discussion, but the concepts can easily be extended to 
more complex physics. 

Eqs. (12. 5p . (12. 6p are constraints which are satisfied at all times provided they are 
satisfied at initial time. When dealing with approximations of the Euler-Maxwell system, 
we will have to check that the proposed schemes are actually consistent with discrete 
versions of these constraints. So, from now on, we will separate these constraints from the 
main system. That they are satisfied will be a property of the proposed approximations. 

The scaling of this problem as well as the material contained in this section is developed 
in detailed in [21] . Here, we briefly summarize the scaling hypotheses. The scaling of the 
hydrodynamic unknowns is the same as in the electrostatic case (see section II. 2p . In 
addition, the charge and current density scales are fixed to po = cuq and jo = eriQUQ. 
The electric field scale is fixed in a consistant way with the electrostatic case, namely 
Eq = Po/ {exoUo). The magnetic field scale is fixed by Bq = Eq/uq. With this choice of 
scales, two dimensionless parameters remain: the scaled Debye length A and the ratio of 
the velocity scale to the speed of light: 

a = —, A = 2 • (2-7) 



c 



We link a and A in such a way that the resulting limit model keeps the largest number 
of terms (the so-called least-degeneracy principle). This principle gives a = A. This 
collection of scaling hypotheses and principles gives rise to the scaled Euler-Maxwell 
model as written above. 
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2.1.2 The A — 7- limit and the Quasi-Neutral Euler-Maxwell system 

We now investigate the formal limit A — )■ of the scaled EM system. 

Proposition 2.2 The formal limit A — ?■ o/ the scaled EM system is the Quasi-Neutral 
Euler-Maxwell ( QN-EM) system: 

V ■ u° = 0, (2.8) 

dtu'^ + V ■ (m° ® M°) = + M° X 5°), (2.9) 

9t5° + VxE° = 0, (2.10) 

-Vxfi° = u°, (2.11) 

V-5° = 0, (2.12) 

n° = 1, (2.13) 



Proof: The formal limit of the EM system is: 

dtn° + V ■ (nV) = 0, (2.14) 

9t(nV) + V ■ (n°M° ® M°) + Vp(n°) = -n\E° + x S°), (2.15) 

dtB^ + V X E° = 0, (2.16) 

-V X 5° = nV. (2.17) 

By assumption, the sequence of initial data satisfies the constraints (12.51) . (12. 6p and con- 
sequently and satisfy (I2.12p . (I2.13P at time t = 0. From the divergence of (I2.17p . 
we deduce (12.81) and that dtn^ = 0. Therefore, (12.131) is satisfied at all times. The other 
equations (EH), (l230ll . ( 1211]) follow immediately ■ 



In this model, the divergence free constraint on is a consequence of (12. lip , while the 
divergence free constraint on 5° is a consequence of (12.101) (and of the divergence free 
initial data). The time evolutions of and 5° are constrained by (12. lip . E'^ is the 
Lagrange multiplier of this constraint. If we resolve this constraint, we find the following 
model: 

Definition 2.3 The reformulated QN-EM model (or RQN-EM model) is: 



V-M° = 0, (2.18) 

dtu'' + V ■ (u° ® u°) = -(E° + u° X E°), (2.19) 

dtB° + V X E° = 0, (2.20) 

V X (V X E°) + E° = -V ■ (m° ® M°) - M° X 5°, (2.21) 

V ■ 5° = 0, (2.22) 
n° = 1, (2.23) 
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Proposition 2.4 The RQN-EM model is equivalent to the QN-EM model provided that 
the Ampere eq. ^2.11\) is satisfied at time t = 0. 



Proof: Let (r2°, n°, 5°) be a solution of the QN-EM model. Taking the curl of fl2:T0|) . 
adding it to fl2.9p and using (12. lip to cancel the time-derivatives leads to f l2.2ip . Con- 
versely, if (n", M*^, is a solution of the RQN-EM model, by proceeding in the same 
way backwards, we easily find that dtiV x 5° + m°) = 0. Therefore, (12.111) is satisfied for 
all times as soon as it is satisfied initially. ■ 

Eq. (12.2 ip is a well-posed elhptic equation for (provided suitable boundary conditions 
are given, such as perfectly conducting or absorbing boundary conditions). In the QN- 
EM model, the hyperbolic character of the Maxwell equations is lost: E^ adjusts to the 
variations of instantaneously. If the initial conditions of the EM model do not satisfy 
(12.111) . an initial layer occurs, during which high frequency oscillations are produced. The 
QN-EM model produces some kind of time averaging of these high frequency oscillations. 

Remark 2.1 If we neglect the inertia of the electrons, which amounts to removing the 
drift term in the momentum equation ^2. 9^ . the QN-EM model reduces to: 

dtB^ + V X E° = 0, 
n° = -V X 
E° + n X B° = 0, 
V ■ 5° = 0, 

which is the so-called Electron- MagnetoHydrodynamics (EMH) system JJJjj. Here, we do 
not make any assumption about the electron time scales, which leads to a slightly more 
complex dynamics. 

In the limit A — t- 0, a change of type of the Maxwell models occurs: it shifts from 
hyperbolic to elliptic. This is the signataure that the EM model is a singularly perturbed 
problem in the limit A — )■ 0. In the process of building an AP scheme, the first step 
is to reformulate the EM model in such a way that this singular perturbation character 
appears more explicitly. This leads us to introduce the following: 



Definition 2.5 The Reformulated Euler-Maxwell (REM) model is: 

dtu^ + V ■ (nV) = 0, (2.24) 

dt{n^u^) + V • {n\^ ® u^) + Vp(n^) = -n^{E^ + x 5^), (2.25) 

dtB^ + V X = 0, (2.26) 
X^d^E^ + V X (V X E^) + n^E^ = -V ■ (n V ® u^) - Vp(n^) - n V x B^, (2.27) 

V ■ 5^ = 0, (2.28) 

A^V-E^ = p^ := 1 -n^, (2.29) 
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Proposition 2.6 The REM model is equivalent to the EM model provided that the Am- 
pere eq. ^2.4^ is satisfied at time t = 0. The formal limit A — )■ 0/ the REM model is the 
RQN-EM model. 

Proof: We proceed like for the proof of proposition 12.41 We take tlie curl of (12. 3p . add it 
to (12. 2p . and use (12. 4p to eliminate the time derivatives of nu and B. This leads to (I2.27p . 
Conversely, proceeding backwards leads to dt{\'^dtE^ — V x — n^u^) = 0. Therefore, 
(12. 4 p is satisfied for all times as soon as it is satisfied initially. The second sentence of the 
proposition is obvious. ■ 

Eq. (I2.27P is a wave equation for E^ with wave-speed A^^. The condition that (12. 4p must 
be satisfied initially provides the Cauchy datum on dfE requested by this time second 
order problem. The use of the REM model preferably to the EM model, in conjuction 
with an implicit time discretization of (I2.27p . is the key for the build-up of an AP scheme 
for the EM model in the quasi-neutral limit A — )■ 0, as we will see in the next sections. 

2.2 Time Semi-Discretization, AP property and linearized sta- 
bility 

2.2.1 Time-Semi-discretization 

The classical scheme for Euler-Maxwell equations uses a semi-implicit discretization of the 
Maxwell equations otherwise the scheme for the Maxwell part is inconditionally unsta- 
ble. The stability requirement of S. Fabre [3l] extended to the electromagnetic case also 
requests the Lorentz force in the momentum equation to be implicit. As shown below, 
this classical scheme requires a CFL condition of the type 6 < 0(A). The AP scheme 
will require two additional levels of implicitness: the first one is an implicit evaluation of 
the mass fiux, like in the electrostatic case (see section [TT^ . The second one is a totally 
implicit scheme for the Maxwell part. The classical and AP schemes can be put in a 
unified framework in the definition below 

Definition 2.7 Introducing the following discretizations: 

^-l(^A,m+l _ ^X,rn^ ^ y . ^^x,m+a^X,m+a^ ^ ^2.30) 
^-l(^^A,m+l^A,rn+l _ ^A,m^A,m^) ^ y . (^A,m^A,m ^ ^\,m^ ^ Vpiu^'"^) = 

_^^A,m,+l-a^A,m+l ^ ^A,m,^A,rn ^ ^\,m-^^ ^2.31) 

^-l(^A,m+l _ ^A,m^ ^ y ^ ^A,m+6 _ ^2.32) 

y2^^l^^\,m+l _ ^X,m.^ _ y ^ ^A,m+c ^ ^\m+a^\,m+a ^ (^2.33) 

with a, b and c taking the values or 1, we define: 

(i) The classical time semi-discrete scheme: {a,b,c) = (0,1,0) or (0,0,1). We will 
consider the case {a,b,c) = (0,0,1) (i.e. the implicitness is in the Ampere eq.) to 
be specific. 
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(a) The AP Semi-Discrete Euler- Maxwell scheme (SD-EM): {a,b,c) = (1, 1, 1). 



For both schemes, the discrete Gauss equation and the divergence free constraint on B 
are satisfied: 

Proposition 2.8 For both the classical and SD-EM schemes, we have: 

V ■ 5^'™ = 0, (2.34) 
A^V ■ E^'™ = 1 - n^'", (2.35) 

for all m G N, assuming that these equations are satisfied initially (i.e. assuming that the 
initial data satisfy 1^2. 34\) , ^2. 35\) with m = 0). 



Proof: Taking the curl of ^M), we find that V ■ B^-"" = V ■ for all m. Eq. 
(12.341) follows from the assumption on the initial condition. To prove (12.351) . we take 
the divergence of (l233l) and eliminate V ■ (n^'"^+''u^'"'+^) thanks to (KM . We find that 
A^V • E^'"" - 1 + n^'™ = A^V ■ E^'° - 1 + Again, follows from the assumption 

on the initial condition. ■ 

We note that the mass fiux in the mass conservation equation and the current in the 
Ampere equation must have the same degree of implicitness in order to guarantee the 
consistency with the Gauss equation. For the SD-EM scheme, it is convenient to use an 
explicit evaluation of the density in the Lorentz force (I2.3ip . because this reduces the 
complexity of the inversion of the implicit scheme. This choice does not restrict the AP- 
character of the scheme nor does it change its linearized stability properties. The classical 
scheme cannot be AP because, when taking the limit A — )■ 0, it does not provide a valid 
recursion for the computation of the variables at time m + 1, knowing their values at time 
m. By contrast, the SD-EM scheme does provide a valid recursion. 

2.2.2 Reformulation and AP property 

The SD-EM scheme can be reformulated into the following scheme: 

Definition 2.9 The Reformulated Semi-Discrete Euler-Maxwell scheme (RSD-EM 
scheme) is: 

^-l^^A,m+l _ ^X,m^ ^ y . ^^x,m+l^X,rn+l^ ^ ^2.36) 

_(^A,m^A,m+l ^ ^X,m^X,m ^ ^X,m^ ^ ^^.37) 
^-l(^A,m+l _ ^X,ra^ ^ y ^ ^X,ra+1 ^ ^2.38) 

(A^ + 6^n^'"')E^'"'+^ + x (V x E^'"+^) = 

-6\V ■ (n^'"*M^'™ ® u^'"") + Vp(n^'™) + n^'"'u^'"' x 5^'™), (2.39) 
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Proposition 2.10 The SD-EM and RSD-EM schemes are equivalent. 



Proof: From the SD-EM scheme (see definition 12.71 with (a, 6, c) = (1,1,1)), we insert 
( I2.32P and f l2.3ip into fl2.33p and get f l2.39p . Proceeding backwards, we get the equivalence 
of the two schemes. ■ 

Now, we investigate the limit A — )■ 0. We first state: 

Proposition 2.11 The formal limit X ^ of the RSD-EM scheme (with 'well-prepared' 

initial data satisfying the constraints 11^2. 34\ ), i\2. 35\) at m = for all X, and such that 

their limits satisfy V x B^'^ + = 0) is the following Reformulated Semi-Discrete 
Quasi-neutral Euler-Maxwell scheme (RSD-QN scheme): 



V ■ (u°'™+i) = 0, (2.40) 

_(^0,m+l ^ ^0,m ^ ^0,m^^ ^2.41) 

^-l(^0,m+l _ ^O.m^) y ^ ^0,m,+l _ ^2.42) 

^o,m+i ^ y X (V X E°''"+^) = -(V ■ (u"'"" ® u"'™) + u^'"" X 5°'™), (2.43) 

V ■ 5°'™+^ = 0, (2.44) 
n°'™ = 1. (2.45) 



Proof: The formal A — of the RSD-EM scheme is: 

^-l(^0,m+l _ ^0,m^) ^ y . ^^0,m+1^0,m+l^) _ ^2.46) 

_(^^0,m^O,m+l ^ ^0,m^O,m ^ ^O.m^,^ (^2.47) 
^-l(^^0,m+l _ ^0,m^ ^ y ^ ^0,m+l ^ ^2.48) 

-6{V ■ (n°''^M°''" ® M°'™) + Vp(n°''^) + n°'"*M°'™ x 5°'™). (2.49) 

Now, taking the divergence of fl2.47p and substracting times the divergence of (12.490 
leads to V ■ (^n^'^+^u^'^+^'j = Q and to n'^'"^+'^ — n^'"^ = 0. Using that the sequence 
of initial data satisfies (I2.35P for all A, we deduce that n^'™ satisfies the initial condition 
^0,0 _ ]^ from which n"'*" = 1 for all m follows. Inserting this into fl2.46p through (12.490 
leads to fl^lO]) , dMU), (1^:^21) and to: 

-5(V ■ (m°'"' ® M°'™) + M^'"* X 5°'™). (2.50) 

But, adding (I2.4ip to the curl of (12.421) and substracting times (12.501) leads to V x 
j^o,m _|_ ^o,m _ g f^j, m G N*. By assumption, the initial condition also satisfies this 
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condition, with m = 0. Inserting it into f l2.50p leads to (12.431) . The divergence free 
constraint on i^o-^+i follows from taking the divergence of f l2.42p and the divergence free 
assumption on the initial data. ■ 



Proposition 2.12 The RSD-QN scheme is consistent with the RQN-EM model. 

Proof: Obvious by comparing the RSD-QN scheme and the RQN-EM model. ■ 

Propositions 12.11] and I2TT2] show that, in the limit A — ?■ 0, the SD-EM scheme is consistent 
with the Quasi-Neutral Euler-Maxwell model. In other words, the SD-EM scheme is AP, 
provided that it is Asmptotically stable. In the section below, we show that it is so, at 
least in the linearized setting. 



2.2.3 Stability analysis 

We will not provide the details of the stability analysis, because it follows the same 
principles as exposed in section 11.41 The details can be found in [21] . Again, to mimic 
the effect of a space decentering in the hydrodynamic part of the model, we consider 
a viscous Euler-Maxwell model, with viscosity f3 satisfying f ll.53p . We also consider a 
linearized model about the state = 1, m'^ = 0, = 0, = and proceed to an 
stability analysis in Fourier space. Again, to mimic the effect of the space discretization, 
we restrict ourselves to Fourier modes satisfying |^| < n/h, where h is underlying spatial 
discretization. Again, the Fourier transform of the solution to the linearized viscous EM 
system can be written in the form fll.6ip . We use definition 11.121 for stability and we say 
that the scheme is Asymptotically Stable if it is stable under a condition on the time-step 
S independent of A when A — )■ 0. In [21], we prove: 

Proposition 2.13 (i) The Classical Schemes is not Asymptotically Stable. 

(a) The SD-EM scheme is Asymptotically Stable, i.e. it is stable under the CFL condition 

S <rh where T is a constant independent when A for A — t- 0. 



2.3 Spatial discretization: enforcing the Gauss law 
2.3.1 Classical fully-discrete scheme 

We present the spatial discretization in the one- dimensional framework. In this frame- 
work, we keep the x and y components of the velocity and of the electric field, and the z 
component of the magnetic field. 

The discretization follows the same general principles as for the Euler-Poisson case 
(see section [T3!) . The discrete unknowns are n|™, u^l"^, Uyl"^, Ex\^^^^2y ^z\'k+i/2^ ^y\T- 
We denote by /n|fc!|_i/2) /Mxlfc+i/2' /«ylfc+i/2 exphcit hydrodynamic fluxes for the mass 
and X and y-components of the momentum conservation equations respectively, while the 
implicit mass flux of the AP scheme will be denoted by /n|^i/2 
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Definition 2.14 The classical fully discrete scheme is: 

5~\n\r' - n\T) + - /niri/2) = 0, (2.51) 

= -n\r'E.\r'-inny)\^B^\T, (2-52) 

= -n|r'i?,ir' + (™.0ir5.|r, (2.53) 

^"'(5^ir+V2 - ^^ir+1/2) + h~\Ey\t_,, - = 0, (2.54) 

^(-^x|™^^l/2 ~ -^xlfcVl/2) = /ri|fc+l/25 (2.55) 

\'6-\Ey\-^' - EyY^) + /i-^(5.|r+Y/2 - fi.iri/2) = (2.56) 



where 



Bz\T — iji.Bz\k'+i/2 + -^zir-1/2)' Ex\^^^ — -{,Ex\'^;^y^ + Ex\^^y^). (2.57) 



r/ie numerical fluxes /«Jr+i/2 ^'^^ /%lfcVi/2 ^'^^ i'^^^'^ ^V- 



fn 1 


m 

fc+1/2 




|m 

lfc+1/2 


/"a 


1 m 

lfc+1/2 




\m 
lfc+1 



(nt.2+p(n))|^ I + I (n«2+p(n))|-, | + 



+^^l\l/2 I (nn,.)|r-('^^x)lr+i I ^ (2.58) 
{nuy)\f - {nuy)\^^. 




where /i^;^/2 ^'^ viscosity matrix. 



Again, we do not make any specific choice of a viscosity matrix and refer to the cited 
hterature for typical expressions of the viscosity matrix. Like in the semi-discrete scheme, 
the classical scheme cannot be AP, because, in the limit A — )■ 0, it does not provide a 
valid recursion to compute the unknowns at the time-step m + 1. We now describe the 
AP fully discrete scheme. 
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2.3.2 AP fully-discrete scheme 

Definition 2.15 The AP Fully- Discrete Euler-Maxwell scheme (FD-EM scheme) is: 

5^\n\r' -n\^)+ /.-^(/„irj/2 - /nir-v2) = o, (2.59) 

= -n\TE.\r'-inUy)\TB.\T, (2.60) 



= -n\^Ey\^^' + {nu,)\^B^\^, (2.61) 

5""^(52ir+l/2 - Bz\k+i/2) + ^""^(-^yir+l"^ - -^■yl™'^^) = 0, (2.62) 

^{Ej:^^;^^!^ - -ExlfcVi/2) = /nir+i/2' (2.63) 

A^ri(E,ir' - Ey\-) + /."^(i?.irJ/2 - BAt%) = {nuy)\r\ (2.64) 



where 



fj \m ^ t D I'm I D I™ ^ Z? im+l ^ f rp im+l i rp im+l \ f-p-N 



The numerical fluxes fuJ]^^i/2 ^'^^ /Mj;|fcVi/2 ^'^^ S'^'ven / (^. Jgj) . while the implicit flux 

k+l/2 



fn[^t^n ^5 9iven by: 



f \m+l _ f \m _ („\m rp |m+l 

J"lfc+l/2 ~ Jn\k+l/2 2 '■Ifc+l '''Ifc ^ -^2;|fc+i/2 

-^(/.jr+3/2 - /.jr-1/2) - + i?.ir+i/2. (2.66) 



Like in the Euler-Poisson case, fl2.66p is obtained from impliciting the central discretization 
part of the classical flux f l2.58p . while keeping the numerical viscosity term explicit. Then, 
using the momentum balance equation fl2.60p and the same kind of lumping of the electric 
field average as for f ll.73p . it is easy to derive f l2.66p . The details are left to the reader. 

The current {nu)x in the x-component of the Ampere equation fl2.63p is evaluated by 
using the mass flux fn\k+i/2- level of the continuous problem, these two quantities 

are identical. Therefore, this approximation is consistent. However, using the mass flux 
rather than the current allows us to guarantee a perfect consistency with the Gauss 
equation, as shown below: 

Proposition 2.16 The Gauss equation 

\'h~\EX^,/, - EX-i/2) = 1 - (2.67) 

is satisfied for all m. G N* provided that it is satisfied by the initial condition (i.e. with 
m = 0). 
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Proof: Taking the difference of f l2.63p evaluated at Xk+1/2 and Xk-1/2 and using (12. 59 p . 
we easily check that: 

X'h~\E,\-'^^^^ - + = \'h~\E^\^^,/, - E^\^_,„) + n|r. (2.68) 

Then, proposition 12. 161 follows. ■ 

We note that this proof would apply to the classical scheme as well. In the ?/-component 
of the Ampere equation (12.631) . the current is evaluated using the usual approximation 
{nuy)\^^°' because, in a one-dimensional problem, the ^/-component of the mass flux is 
independent of x and does not enter the mass balance. In a 2 or 3-dimensional problem, 
one should evaluate all components of the current using the corresponding components of 
the mass flux, to ensure consistency with the Gauss equation. 

We now reformulate the FD-EM scheme in order to study the limit A — !■ 0. 



2.3.3 Reformulation and AP property 

Definition 2.17 The Reformulated Fully-Discrete Euler-Maxwell scheme (RFD-EM 
scheme) is 

ri(n|r^ - nlD + - /nirv2) = 0> (2-69) 

ri((nn,.)ir' - {nu^rn + h^\UT^i/2 " = 

= -n\^E^\r'-{nUy)\^B,\T, (2.70) 
S-\{nUy)\r' - {nny)\T) + /^-'(/.jr+1/2 " fuS-1/2) = 

= -n\^Ey\^^' + inu.,)\^B^\T, (2.71) 
S-'iB.\T:i/2 - B.\T+i/2) + h-\Ey\^^,' - Ey\l^^') = 0, (2.72) 



(A' + -wi^lT+i + n\T)) i^.ir+Y/2 = A'^.ir+1/2 + Sfn 
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2 I'^lfc+l "Ifc )J ^^\k+l/2 ~ ^ ^x\k+l/2 T" "./nlfc+1/2 

2 (/«^lfc+3/2 ~ fuAk'^1/2) ^(("-^j/)!" + ('^^y)lfc+i) -^2lfc+i/2' (2.73) 
(A^ + 5^n|^)E,ir' - ^'h-\Ey\^ti - 2^.ir' + Ey\t^l) = ^'Ey\T + 6{nUy)\T 

-^h^'mt^,,, - 5,1-1/2) - ^'h~\us^,/, - /„jr-i/2) + nnu.)\k bm\ (2.74) 

with Ex\^ and -Bzl™ given by Ii2.65\) . 

Proposition 2.18 The FD-EM and RFD-EM schemes are equivalent. 

Proof: We first consider E"^^^: inserting (12.621) and (I2.6ip into (12.641) to eliminate 
-^2|fc+i/2' ('^Wy)ir'^^ respectively, we find that (I2.64p is equivalent to (I2.74p . 

We now examine E"^^^ . To get (I2.73p . we simply insert the expression (I2.66P of the 
mass flux into ( 12.63p . ■ 
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Eq. f l2.74p is a well-posed discrete elliptic equation provided that suitable boundary 
conditions are defined. It allows to compute Ey\^^^ from known values. Eq. fl2.73p 
provides a direct explicit evaluation of E^l^^^^. 
We now investigate the A — limit. We have: 

Proposition 2.19 (i) The formal A — )■ limit of the RFD-EM scheme (with 'well- 
prepared' initial data satisfying the discrete Gauss eq. {2.67^ at m = for all A and 
such that their limits satisfy Uy\l — h~^{B ^l^j^^^^ ~ -^2|fe-i/2) = Oy) «s 

= 1, (2.75) 

1/ im+l |m.\ I U — ^ff I'm r \m \ 771 im.+l „ I'm 75 im fry ^(;\ 

^ {^x\k ~ ^x\k ) -r yJu^\k+l/2 ~ Jua:\k-l/2) — ~'^x\k ~^y\k^z\ky {■^■10) 

6~\uy\r' - uy\T) + - fujT-1/2) = -Ey\r' + (2.77) 

^'\B,\T:i/2 - B-AT-,1/2) + h-\Eyra,' - Ey\^+') = 0, (2.78) 

771 im+l ( f \^ f I™ —{^1 I™ _L 1/ I™ R I™ _L 

=^lfc+l/2 ~ 2 ^■'^^\k+3/2 JUx\k-l/2) 2^ "T "ylfc+lJ -'^2|A:+1/2 "T 



+'5"Vnir+i/2, (2.79) 

771 im+l A^^/pi im+l _ op im+l , rp \m+l\ 

^y\k \^y\k+i ^^y\k ' -^ylk-i ) ~ 

= —h ^{fuy\'k'+l/2 ~ fuy\'k-l/2) + ^x\T Bz\T- (2.80) 

with Ex\k ^''^'^ ^z\k fl'^'^^en hy Ii2.65\) and with the numerical fluxes: 

fn\k+l/2 \ 1 / M^;!^ + \ u"*^^ ,„ / 

(^M^ We denote hy: 



/^r+i/2 I UxYk — Ux^^^i I • (2.82) 





m 

fc+1/2 




hn 

1 fc+1/2 




Im 

lfc+1/2 



Then, for m > 1, ([^. yg[) zs equivalent to: 



m+l m m ^ ,n m 

-^a:^lfc+l/2 ^ ^(/«xlfc+3/2 ^ 1^-1/2) ~ ■^(^?/ir + %lfe+l) 



2 

1 



fc+1/2 



/ 771 im 9771 |m \ T? \m 

~^l-^x|/c+3/2 ~ ^-'^2;|fc+l/2 + ^x\-k- 



1/2) 



^ („. \m~l(-r> im— 1 p im-l \ l'"~l('R D im— 1 \\ 

^l"ylfc+i l-Dzlfc+3/2 -°2lfc+i/2J "ylfc l-°2lfc+l/2 -°^lfc-l/2J>' 

+ — (-Dnir+1/2 - -Dn.ir+1/2)- (2-83) 

for all m > 1. Then, this scheme ( the Reformulated Fully-Discrete Quasi-Neutral Euler- 
Maxwell (RFD-QN) scheme) is consistent with the RQN-EM model. 
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Proof: (i) Taking the limit A in the RFD-EM scheme only modifies (12T3D . ((221) 
into 



S 

2^"'\k+l "I" "'Ik ) -^a;|fc+l/2 ~ Jn\k+l/2 



- UT-1/2) - + i?.ir+i/2, (2.84) 

-/."^(i?,|r+i/2 - 5.ir-i/2) - 5/^^'(/.jr+i/2 - Ajri/2) + 6{nu,)\TB^\T. (2.85) 

Comparing (12.841) and fl2.66p shows that fn\^^i/2 = 0, and, inserting it into (12. 69 p . that 
(I2.75P is satisfied for all m > 1. With the assumption that the initial condition satisfies 
( I2.67P for all A, we deduce that n\'^ = 1 for all integer m > 0. Inserting this into fl2.70p . 
fl2TB . (12T2|) . (EUD and ( 1238|) . we are led to (l2T6|) . ( K77\f . fl2T8D . (I279|) and (ICT]) . 
Concerning fl2.85p . at this point, we get: 

rp im+l h~'^('P _ 9 Z? I p |ni+l\ 

-^ylfc \^y\k+i ^^y\k ~^ ^y\k~i ) ~ 

U — iff ym f \rn ^i,, ['"Rl™ 

— \Juy\k+l/2 Juy\k~l/2J ^ "-xlk -^zlk 

+r^K|r - h^'mk_,,/, - 5.1^-1/2)). (2.86) 

However, taking the differences of (12.780 for A; + 1/2 and k — 1/2, multiplying the result 
by and subtracting the result to (12.770 leads to 

UylT^^ - ^"^(^^1^+1/2 - ^^ir^i/2) = 0' (2.87) 

for all m > 0. But with the well-prepared initial data assumption, (I2.87P is valid with 
m + 1 repaced by m, for all m > 0. Inserting this identity into (I2.86P leads to (I2.80p . This 
concludes point (i) of the proposition. 

(ii) All equations in the resulting scheme are consistent with the RQN-EM model, 
except (I2.79p . where an 0{6~^) terms appear: <^~Vn|fc!|_i/2- ^^^^ P^^^ proof, we 

show that fn\'k+i/2 can be substituted with the last three fines of (I2.83p . Indeed, taking 
half the sum of (12.761) for k and k + 1, we can write: 

1 1 

r-i/zf im+i _ _n ( f \^ _r) I"* "i^i 

\.\Jn\k+l/2 2 ''=+1/2'' Un|fc+l/2 2 "'^+1/2'''' ~ 

f r \m r \m \ 

~ 2 y-'^x'k+3/2 JUx\k-l/2) 

\( fP [m+l I 9 p \m+l I 771 |m+l \ 

^\^^\k+3/2 '^^^\k+l/2 ^^\k-l/2! 

~^('"ylfc+i(-Sz|fc+3/2 + -^21^+1/2) + '^y\k'{Bz\k'+i/2 + -^21^-1/2)). (2.88) 



37 



The, subtracting fl2T9D to (12:881) leads to 



e— 1 /■ im+l — ( \m+l op |m+l i p N 

'^\"'y\k+l\-'^z\k+3/2 ^z\k+l/2) "-ylk \^z\k+l/2 ^z\k~l/2)) 
+ ^Pn|r4"v2-^n|r+l/2)- (2-89) 

Substituting 5~^/n|fc!|_i/2 (I2.79P by its expression deduced from fl2.89p (with m replaced 
by m — 1) leads to fl2.83p for all m > 1. 

Now, we see that the second and third lines of fl2.83p are at leading order, equal to 
-(/iV4) ((92E,.)lr+i/2 + {d^{uyd^B Mk7i/2) and tend to zero with h. 

Then, the fourth line can be expanded using Taylor's expansion and is of the order of 
{dtDn){xk+i/2-,t^~^^'^) (supposing that the artificial viscosity term can be interpolated by 
a smooth function Dn{x,t)). But the artificial viscosity is 0{h) (see (I2.82P ). Therefore, 
this term is 0{h) and tends to zero with h as well. 

Therefore, fl2.83p is consistent with f l2.2ip (having in mind that, in this 1-D geometry, 
the X component of V x (V x i?) is identically zero. Then, the RFD-QN scheme is clearly 
consistent with the RQN-EM model. This concludes the second point of the proposition. 



This proposition shows that the RFD-EM scheme is AP. The assumption of well-prepared 
initial conditions can be removed. In this case, the scheme takes its form as stated in 
the proposition for m > 1. We note that, in the Euler-Maxwell case, there is no need for 
an inverse CFL condition when the viscosity terms -0^1 ^^^^2 not identically zero (see 
discussion at the end of section [1.5.41 for comparison with the Euler-Poisson case). 



2.4 Euler-Maxwell problem: conclusion 

In this section, following the same strategy as for the Euler-Poisson problem, we have 
provided an Asymptotic-Preserving (AP) scheme for the one-fluid Euler-Maxwell problem 
in the quasineutral limit. However, the quasi-neutral limit of the Euler-Maxwell model 
is more complicated than that of the Euler-Poisson problem and gives rise to infinite 
propagation speed of electromagnetic waves. Again, following the previous strategy, we 
provide a reformulation of the Euler-Maxwell model which is not singular when A — 0. 
The AP scheme is therefore designed to mimic this reformulated Euler-Maxwell system. 
Again, it can be based on any classical shock capturing scheme and consists in perturbing 
the standard hydrodynamic fluxes by corrective terms which are of the order of the time 
step or the mesh step. However, additionally, a fully implicit treatment of the Maxwell 
equations must also be implemented. A linearized stability analysis has been reviewed. It 
confirms that the stability condition of the AP-scheme is independent of A when A — )■ 0. 
Again, non-linear stability results are still open. 
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3 Extensions 



Pull Euler, Navier-Stokes, etc. The AP-methodology extends straightforwardly when 
the isothermal or isentropic gas dynamics model is replaced by a full Euler system includ- 
ing an energy equation. This extension is straightforward in both the Euler-Poisson and 
Euler-Maxwell cases. Similarly, viscosity or heat conductivity terms can be considered 
without altering the principles of the method. 

Two fluids or more. The case of two-fluid models, where each of the ion and electron 
species is modeled by its own Euler system of equations, and are coupled to the electric 
potential or electro- magnetic field by the electrical sources (charges and currents), has 
been considered. References [T3j and [H] for the Euler-Poisson case, and [21] for the 
Euler-Maxwell case show that the approach extends easily to this case (and would also 
apply to the multiple ion species case as well). The numerical results have been obtained 
in this case, with the physical electron to ion mass ratio. We refer to the refereces for 
more details. 

Euler-Poisson-Boltzmann. A commonly used approximation in plasma physics is to 
suppose that the electrons follow a Boltzmann law. The Boltzmann law provides a linear 
relationship between the electric potential and the logarithm of the electron density (or 
chemical potential). It is obtained from the electron momentum conservation equation in 
the limit of vanishing inertia. Then, the Euler-Poisson-Boltzmann (EPB) system consists 
of a pressureless gas dynamics model for the ions with an electrical forcing. The electric 
potential is obtained by solving the Poisson equation where the electric charge takes 
into account the ion and electron densities, the latter through the Boltzmann law. The 
resulting Poisson equation is nonlinear. In the quasineutral limit, the EPB model reduces 
to the compressible gas dynamics equations, the electrical force acting as a pressure term 
for the ions. The AP methodology has been applied to the EPB model in [22] . 
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Part 2 
Large magnetic fields 
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4 The isentropic Euler-Lorentz model 



4.1 Introduction 

This section and the following ones are concerned with the numerical approximation 
of the Euler equations for charged particles subject to the Lorentz force (the 'Euler- 
Lorentz' system), when the magnetic field is large, or equivalently, when the parameter 
r representing the reciprocal of the non-dimensional cyclotron frequency tends to zero. 
In this regime, the so-called drift-fluid (or gyro-fluid) approximation is obtained. In 
this limit, the parallel motion relative to the magnetic field direction splits from the 
perpendicular motion. The latter is given by an algebraic relation which describes the 
various drifts of the fluid across the magnetic field lines. The parallel motion is given 
implicitly by the constraint of zero total force along the magnetic field lines. In these 
sections, we construct Asymptotic-Preserving (AP) schemes which give rise to both a 
consistent approximation of the Euler-Lorentz model when r is finite and a consistent 
approximation of the drift limit when r — )■ 0. Above all, they do not require any constraint 
on the space and time steps related to the small value of r. 

4.2 Setting of the problem 

The Isentropic Euler-Lorentz (lEL) model consists of the system of isentropic Euler equa- 
tions subject to the Lorentz force. In this study, the electro-magnetic field is supposed 
given with an arbitrary spatio-temporal dependence. Of course, we have in mind that 
it satisfies the Maxwell system or any system derived from it, but the only information 
that we shall use from it is that the magnetic field is divergence-free. For the same rea- 
son, a single plasma species is considered (the ions) but again, all considerations below 
would extend to the electron species and any kind of coupling between these two species 
(either through Poisson's equations or through Maxwell's equations, or again, through 
quasi- neutrality) . 

In this framework, the lEL model is written: 



where, n{x,t) > and u{x,t) G M.^ stand for the ion density and ion velocity and de- 
pend on the space-variable x G and on the time t > 0. The electro-magnetic field 
{E{x,t),B{x,t)) G X is supposed given and satisfies V ■ -B = 0. The positive ele- 
mentary charge is denoted by e and the ion mass, by m. The ion pressure p is supposed 
to be a given function of the density (e.g. p{n) = kBTn in the isothermal case, where T 
is the ion temperature and fc^, the Boltzmann constant, or p[n) = CrC , where 7 > 1 and 
C > are given constants, in the isentropic case). 

The following scaling allows us to highlight the large magnetic field regime. We denote 
the scaling units for length, time, velocity, density, pressure, electric field and magnetic 
field by Xq, to? ""O) ''^O; Po? Eq, Bq. As usual, we relate the velocity scale to the time 



dtn + V ■ (nu) = 0, 

m{dt{nu) + V ■ {nu ® u)) + Vj9(n) = en{E + u x B) 



(4.1) 
(4.2) 
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and space scales by xq = uoto- We also relate the electric and magnetic field scales by 
£^0 = uqBq. This relation indicates that the typical electric field in the plasma is of the 
same order as the electric field induced by the motion of the plasma accross the magnetic 
field lines. We introduce the ion sound speed = (po/ (mno))^^^. The ion gyro- frequency, 
i.e. the angular velocity of the gyration motion about the magnetic field lines is given by 
(jJq = eBo/m = eEo/{muo)- Two dimensionless parameters appear: the Mach number M 
and the scaled gyro-period r given by 

M=^, r = J- = ^, (4.3) 

This leads to the following scaled lEL model: 
dtn + V ■ (nu) = 0, 

dt(nu) + V ■ (nu ® u) + -T7^'^p{n) = -n{E + u x B). 

T 

Now, there are two interesting scales: 

1. Case: = r -C 1. In this case, the pressure force is of the same order of magnitude 
as the Lorentz force, but much larger than the inertia terms. Without the Lorentz 
force, this scaling corresponds to the low Mach-number regime. In the presence of 
a Lorentz force, the low Mach-number regime still applies to the parallel motion, in 
a modified form, as we will see below. 

2. Case: = 0(1), r ^ 1. In this case, the Lorentz force is smaller than the Lorentz 
force. This scaling requires that E\\ = 0{t)E± and the parallel dynamics remains 
that of a compressible fiuid. 

From the viewpoint of AP schemes, the second case is simpler to treat than the first one 
and, for this reason, we will focus on the first case. The developed schemes will obviously 
be suitable for the second well. Therefore, from now on, we assume: 

= r, (4.4) 

which leads to the final scaled form of the lEL model: 

Definition 4.1 The scaled Incompressible Euler-Lorentz model (lEL) is: 

dtn'' + V ■ irfu^) = 0, (4.5) 
r {dt{n^u^) + V ■ {n^u^ ® u^)} + Vp(n^) = n''{E + u' xB). (4.6) 



The following notations will be useful: the direction of the magnetic field is denoted by 
h = B/\B\ wherever B 0. Any vector quantity v can be split into its parallel and 
perpendicular parts as follows: 

V = v\\b + v± , v\\ = V ■ b , v± = b X {v X b) = (Id — b ^ b)v , (4.7) 
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where the matrix (Id — 6® 6) is nothing but the projection matrix onto the perpendicular 
plane to b. Next, we introduce the parallel and perpendicular gradients of a scalar function 
(p by 

V||0 = (V0)|| =6- V0, Vx0= (V0)^ = 6x (V0X 6) = (Id-6®6)0. (4.8) 
Similarly, the parallel and perpendicular divergence of a vector field v are defined by: 

V|| ■ w = V ■ (i7||6), V±-v = V-v±, V ■ V = V\\ ■ V + V±- V. (4.9) 
We also note that, since \b\ = 1, any derivative of b is orthogonal to b, i.e. 

b.dtb = 0, 6- = 0, VA; G {1,2,3}. (4.10) 

5 The Drift-fluid Limit r ^ 

In this section, we investigate various formulations of the model obtained by letting r — )■ 
in the lEL model (03]), fOj) . 

5.1 The Isentropic Drift-Fluid (IDF) model; a first reformula- 
tion 

The formal limit r — )■ in the lEL model (14. Sp . (14. 6 p leads to the so-called Isentropic 
Drift-Fluid (IDF) model: 

Definition 5.1 The Isentropic Drift-Fluid model (IDF) is: 

+ V • (n°n°) = 0, (5.1) 
Vp(n°) = n°(E + M° X 5). (5.2) 



We obviously have: 

Proposition 5.2 The formal limit of the lEL model is the IDF model. 
We now propose a first reformulation of the IDF model: 

Proposition 5.3 A solution of the IDF model is a solution to the following reformulation 
(first Refomulated IDF model or (RIDF-1) model): 

dtn^ + V ■ (nV) = 0, (5.3) 

= 4^ X (Vp(n°) - n°E) , (5.4) 
B 

-V|| (p'(n°)V|| ■ (n°nj)) = dt{n^E\\) - dtb ■ Vp{n°) + Vy (p'(n°)Vx ■ (^^°0) • (5.5) 

The converse is true provided that the initial data satisfy the constraint (V||p(r;,°) — 
nOE||)|i=o = 0. 
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Proof: We split (15. 2 p into its parallel and transverse components with respect to the 
magnetic field direction b. First, concerning the transverse component, we take the vector 
product of fl5.2l) with b. The resulting equation can easily be resolved for m*[ in the form 

Taking now the scalar product of fl5.2p with b, we find: 

V||p(n°) - = . (5.6) 

This equation does not explicitly contain u^^ but defines a constraint which indirectly 
determines it. The resolution of this constraint leads to the elliptic equation (15. 5p . To 
show it, we first, multiply (15. ip by p'{n^) and we get: 

atp(n°) + p'(n°) (Vx ■ (n^O + V|| ■ (n°njj)) = . (5.7) 

Applying Vy to (15. 7p . noting that [Vy, dt] = —dtb-V (where [■, •] denotes the commutator) 
and inserting (15. 6p leads to (15. 5p . Reciprocally, it is straightforward to see that system 
(15. 3p . (15. 4p . (15. 5p . implies system (15. ip . (15. 2 p provided that eq. (15.60 is satisfied at time 
t = 0. This ends the proof of proposition 15.31 ■ 

After dividing by n", we find that the first term at the right-hand side of (15. 4p is the 
diamagnetic drift velocity while the second one is the E x B drift velocity. Eq. (15. 5p 
is a well-posed one- dimensional elliptic equation for posed along the magnetic field 
lines, provided that adequate boundary conditions are given. The determination of the 
boundary conditions for (15. 5p will be discussed in more details below. 



5.2 Analogy with the low Mach-number limit of compressible 
fluids 

In this section, we depart from the drift-fiuid limit and consider the low Mach-number 
limit of ordinary isentropic compressible fiuids. We show that the procedure which leads 
to (15. 5p is specific to the Euler-Lorentz problem and cannot be reproduced in the case of 
ordinary fiuids. In the next section, we will use the analogy with the low Mach-number 
limit of ordinary fluids and devise an alternate expression of the limit problem (I5.3p - (l5.5p . 

The scaled Isentropic Compressible Euler (ICE) system with the low Mach-number 
scaling is written as follows 

dtn^ + V ■ (n^u^) = 0, (5.8) 
r {dt{n^u^) + V ■ (n^u^ ® u^)} + Vp(r2^) = 0. (5.9) 

Then, in the limit r — 0, we formally get: 

(9tn° + V ■ (nV) = 0, (5.10) 
Vp(n°) = 0. (5.11) 

We suppose that the boundary conditions are such that (15. lip is well-posed and gives 
p(n°) = Constant. This occurs for instance if the boundary conditions for are mixed 
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Dirichlet or Neumann conditions and if the values along the Dirichlet boundary are uni- 
form. More general conditions are of course possible but will not be detailed here. Here, 
we additionally impose that these conditions lead to the fact that p{n^) and therefore nP 
are independent of both space and time. As a consequence, (15.101) leads to 

V-M° = 0, (5.12) 

However, this is not enough to determine u^. But, dividing (15.91) by r and letting r — 
leads to the existence of a scalar function such that 

dtu° + V ■ (n° ® n°) + Vp^ = 0. (5.13) 

The hydrostatic pressure p^ is the first order (in r) perturbation pressure, i.e. 

pi=(lim^), (5.14) 

and can be determined from the incompressibility constraint (I5.12p . Indeed, taking the 
divergence of (15.131) and using (I5.12p . we find 

Ap^ + : (u° ® n°) = 0, (5.15) 

where denotes the Hessian matrix (matrix of second order derivatives) and : the 
contracted product of tensors. Eq. (15.151) is an elliptic equation which determines p^ 
provided appropriate boundary conditions are given. These conditions can be deduced 
from those for n^. 

As a summary, the low Mach-number limit of the ICE eqs. (15.81) . (15.91) is the Incom- 
pressible Euler (IE) eqs.: 

n° = Constant, (5.16) 
V-n° = 0, (5.17) 
+ V ■ (m° ® M°) + Vp^ = 0. (5.18) 

To identify the limit problem, the strategy of section \5A] would not be not adequate. 
Indeed, if we take the gradient of (I5.10p and use (15. lip , we are led to 

V(V-M°)=0, (5.19) 

which is not a well-posed problem for u^. Indeed, the operator — V(V ■ u) = —An + V x 
(V X u) is not elliptic because of the term V x (V x m), except in dimension 1 where this 
term does not appear. By contrast, in the case of the IDF model, the problem is well- 
posed, thanks to the presence of the Lorentz force. On the one hand, the Lorentz force 
provides the explicit algebraic relation (15. 4 p for u±. On the other hand, u\\ is determined 
by inverting the operator — V||(V|| ■ u\\) (see 15. 5p . which now leads to a well-posed elliptic 
equation because u\\ is a scalar and the equation is posed on a one-dimensional manifold 
(the magnetic field line). 
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5.3 A second formulation of the drift-fluid limit r — > using the 
analogy with the low Mach-number limit 

In the present section, we use the analogy with the low Mach-number limit of ordinary 
fluids developed in the previous section to devise an alternate expression of the limit 
problem fl5.3p -( l53|) . which in turn will be useful for the numerical developments below. 
Specifically, we want to exploit the analogy of the constraint ( 15. 6p with (15.111) . and of 
the continuity eq. (15. ip with (I5.10p (we leave the algebraic relation (15.40 apart because 
it will not play any role in the discussion). In the case of the low Mach-number limit, 
these two equations respectively lead to (I5.16P and to (15.171) . while in the case of the 
drift-fluid limit, they are unchanged, meaning that they cannot be further simplifled. 
However, in the low Mach-number limit, we know that these two relations are not sufficient 
to characterize the limit solution and an additional relation was sought by dividing the 
momentum conservation equation by r and taking the limit. We explore a similar strategy 
here with the parallel component of the momentum equation. 
We first introduce 

T 

= lim . (5.20) 



We have 



:= lim ^^"""^ ^^"""^ = p'in^y (5.21) 

T-s-O T 

Upon dividing (14. 6 p by r and taking the limit r — 0, we deduce that 

b ■ (9j(nV) + V • (riV ® n°)) = -V||(p'(n°)n^) + n^En := F^. (5.22) 

Then, we supplement the IDF model (I5.ip - (l5.2p with the additional eq. (I5.22p and 
introduce the following augmented model: 

Definition 5.4 The 'augmented IDF' model (AIDF model) is : 

V[|p(n°) -r2°E|| = 0, (5.23) 
dtn^ + V ■ (nV) = 0, (5.24) 
6- (at(n°n°) + V-(nV®n°)) = -V \\{p' {n'^)n^) + E\\, (5.25) 

= ^bx (Vp(n°) - n^E) . (5.26) 



We have the obvious: 

Proposition 5.5 Any solution of the AIDF model is a solution of the IDF model. Recip- 
rocally, a solution of the IDF model gives rise to a solution of the AIDF by solving liS. 25\) 
for . 
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The AIDF model allows us to change the viewpoint and instead of seeing fl5.23p as a 
constraint which implicitly determines u^^, we can see it as a constraint which determines 
while u^^ is found by solving the evolution eq. ( ]5.25p . This viewpoint is highlighted in 

Proposition 5.6 Any solution of the AIDF model provides a solution of the following 
model (the second reformulation of the Isentropic Drift-Fluid model or RIDF2 model): 

V||p(n°) - n°E|| = , (5.27) 

-V|| ■ (V||(p'(n°)^^) - ^^^ll) = -d^n^ + V ■ ((6 ® fe)V ■ ® 

-Vx ■ {n^u\dth) - Vji ■ {n\\ ■ dth) - dtV± ■ {n%l). (5.28) 

b ■ {dt(n\^) + V ■ (n%° ® u°)) = -V||(p'(n°)ni) + n^E\\, (5.29) 



o„,o 



1 



n u 



—b X (Vp(n°) - n°E) . (5.30) 



B 



Conversely, any solution of the RIDF-2 model such that the mass conservation eq. ( [5. 24\) 
is satisfied at time t = 0, is a solution of the AIDF model. 

Proof: We proceed like in section 15.21 Taking the time derivative of fl5.24p and using 
f ra . we get: 

+ dtVii ■ (n%J) + dtV^ ■ {n\l) = 0. (5.31) 

We note the identity: 



dtVii ■ {n\?,) = dtV ■ {n\%) = dtV ■ ((6 ® 6)(n°n°)) 



= V ■ ((6 ® b)dt{n\^)) + V ■ {{dtb ®b + b® 9t6)(nV)) 

= V ■ {b{b-dt{n\^))) +V^-{n%ptb) + Vwin\l-dtb). (5.32) 

Multiplying (15.251) by the vector b, inserting it into (15.321) . and using definition (I5.22p of 
we find: 

dtV\i ■ (n%;) = -V ■ {{b ® b)V ■ (nV ® m°)) + Vy ■ F|f + 

+V± ■ {n\ptb) + V|| ■ (n%° ■ dtb). (5.33) 

Inserting (I5.33P into (I5.3ip leads to 

- V ■ ((6 ® 6)V ■ (n°n° ® n°)) + Vy ■ F|f + 



■ (n°Mpt6) + V|| ■ (n°< ■ dtb) + (9tVx ■ (n°0 = 0. (5.34) 



In view of definition (15.220 . (I5.34p leads to (15.280 . Reciprocally, if {nP, u^,n^) is a solution 
of the RIDF-2 model, it is an easy matter to see that it satisfies the AIDF model, provided 
that the mass conservation eq. (15.240 is satisfied at time t = 0. This end the proof of the 
proposition. ■ 
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fl5.28p is a one- dimensional non-linear elliptic equation for posed along the magnetic 
field lines. 

Like in the quasi- neutral limit case, the key for designing AP schemes in the drift-fluid 
limit case is to reformulate the Euler-Lorentz model in the form of a singular perturbation 
of the Drift-Fluid model. In the next section, we provide two different such reformulation, 
which are based on the two reformulations of the Drift-fluid model established above. 

6 Reformulations of the Euler-Lorentz model 

In this section, we construct two equivalent reformulations of the Euler-Lorentz model. 
These reformulations will be the bases for two different AP-schemes. 

6.1 First reformulation of the Euler-Lorentz model: wave equa- 
tion formulation for nu^\ 

The RIDF-1 reformulation of the IDF model has a counterpart at the level of the lEL 
model as shown in the following: 

Proposition 6.1 Any solution of the lEL model is a solution of the following first Re- 
formulated lEL model (RIEL-1): 

dtrf + V ■ {rfu^) = 0, (6.1) 
n^ul = j^b X {VpK) - n^E + r [dt{n^u^) + V ■ {n^u^ ® u^)]} , (6.2) 

\id\ 

rdt (^{dt{n^u^) + V ■ {n^u^ ® n")},,) - V||(p'(r2")V|| ■ K^)) = 

= dt - dtb ■ Vp{n^) + V\\{p'{n^)V^ ■ (n^uD). (6.3) 

Conversely, any solution ot the RIEL-1 model such that Ii6.5\) is satisfied at t = is a 
solution of the lEL model. 

Proof: We apply the same algebraic manipulations to the lEL model fl4.5p . fl4.6p . as we 
did to the IDF one f IS.ip . 05. 2p in section [5?T1 We first split fl4.2p into its parallel and 
perpedicular components. This leads to: 

r {dt{n^u^) + V ■ {n^u^ ® u^)}^ + Vxp(n^) = n^{E^ + ulx B), (6.4) 
r {dt{n^u^) + V ■ (n^u^ ® m^)},, + V||p(n^) = n^^y. (6.5) 

The transverse component eq. (16. 4p can be recast in the form: 

n-'ul = j^b X {Vp(n^) - n'^E + t [dt{n''u^) + V ■ {n^u^ ® u^)]} , (6.6) 
\B\ 

and clearly appears as a singular perturbation of (15. 4p . For the parallel component eq. 
( 16. 5p . we use (15. 7p . which is also valid for finite r. Like in section 15. ![ we apply Vy to 
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i\5.7L commute dtpin^) and Vy and use (16. 5p to eliminate V||p(n'^). This leads to (16. 3p . 
Conversely, it is an easy matter to see that any solution of the RIEL-1 model provided 
that (16.51) is satisfied at t = 0. This ends the proof. ■ 



Eq. (16. 3p is a wave equation for ri^u^^y Indeed, we have, using (I4.10p : 

{dtijiu) + V ■ {nu ® = dt{nu\\) + V ■ {nuu\\) — {dt + u ■ V)b ■ nu±. (6.7) 

Then, (16. 3 P can be rewritten: 

ra«(n"Mp - V||(p'(n")V|| ■ {n^ul)) + rdtV ■ (n^u^ul) = 

= rdtiidt + ■ V)b ■ n^ul) + dt 

-dtb ■ VMn") + V||(p'(n")Vx ■ (n'O). (6.8) 

For small u, the third term of the left-hand side of (16. 8p can be neglected to leading 
order. Then, the principal symbol of the differential operator acting on n'^My is 
'^dtt — j5'(n'')(V||)^, which is a wave operator associated to the acoustic wave speed 
Cs = {p' {n^) / tY ^"^ ■ For the sake of simplicity, let us take an isothermal equation-of-state 
p{n) = nT where T is a fixed temperature. Then, Cs = (T/tY^'^ is the typical velocity of 
acoustic waves in the low Mach-number scaling. 

The requirement that (16. 5p must be satisfied at t = sets up the additional initial 
condition that is needed for the time second order differential equation (16. 3p . 

We note the obvious 

Proposition 6.2 In the limit r — 0, the RIEL-1 formulation of the Euler-Lorentz model 
formally converges to the RIDF-1 formulation of the Drift-Fluid model. 

A first class of AP-schemes for the lEL model will rely on an implicit discretization of the 
wave eq. (16. 3p . 

6.2 Second reformulation of the Euler-Lorentz model: wave 
equation formulation for n 

Now, we show that the RIDF-2 reformulation of the IDF model has also a counterpart at 
the level of the lEL model. By contrast to the previous one, this reformulation involves 
a wave equation for n. To this aim, we apply the methodology of section 15.31 

Proposition 6.3 Any solution of the lEL model is a solution of the following second 
Reformulated lEL model (RIEL-2): 

Tdttrf - V|| ■ (V||p(n^) - E{) = r {V ■ ((6 ® 6)V ■ {n^u" ® u^)) 

-Vx • {rfu\dtb) - V|| ■ {n^u\ ■ dtb) - dtV± ■ {n^uD} , (6.9) 

n^ul = -l-b X {Vp{n^) - n^E + r [dtin^u^) + V ■ {n^u^ ® u^)]} , (6.10) 
{dt{n^u^) + V ■ {n^u^ ® u^)}u + -(V||p(n^) - n^Ey) = 0. (6.11) 
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Conversely, any solution ot the IREL-2 model such that \4-^ ^■^ satisfied at t = is a 
solution of the lEL model. 

Proof: We leave (16 .6^ unchanged. We use f l5.3ip and f l5.32p . which are obviously valid 
for finite r. To eliminate dt{n'^u^) in f l5.32p . we use fl6.5p in the form: 

(6 ® h) {dt{n^u^) + V ■ (n^u^ ® u^)} + -{h ® b){Vp{n^) - n^E) = 0. (6.12) 

r 

This leads to fO]) . 

Reciprocally, it is an easy matter to see that any solution of the RIEL-2 model satisfies 
the original lEL model provided that (14. 5 p is satisfied at t = 0. This ends the proof. ■ 

Eq. (16. 9p is a wave equation for n. If u is small, at the leading order in u, the principal 
symbol of the differential operator acting on n is again rdu — p'(ri'^)(V|j)^, and is again a 
wave operator associated to the acoustic wave speed Cg = {p'{n'^)/TY^'^. 

The requirement that (14. 5 p should be satisfied at t = sets up the additional initial 
condition that is needed for the time second order differential equation (16. 9p . 

Again, we have the obvious: 

Proposition 6.4 In the limit r — 0, the RIEL-2 formulation of the Euler-Lorentz model 
formally converges to the RIDF-2 formulation of the Drift-Fluid model. 

A second class of AP-schemes for the lEL model will rely on an implicit discretization of 
the wave eq. (16. 9p . 

6.3 Boundary conditions 

The question of boundary conditions is complex for several reasons. First, the theory 
of boundary value problems for hyperbolic systems of conservation laws is still in its 
infancy (see e.g. [H [331 160]^ ). Unfortunately, the cases that can be rigorously treated 
seldom apply to practical situations. A practical rule is that the number of boundary 
conditions that can be imposed corresponds to the number of entering characteristics. 
However, this number depends on the solution itself. In practice, the state variables 
are supposed known outside the domain and a boundary Riemann problem is solved 
between the prescribed outer values and the current inner values of the state variables. 
This permits the computation of the entering fluxes and the advancement of the solution. 
Neumann boundary conditions can be prescribed by supposing that the value of the 
corresponding state variable in the outer cell is equal to the value in the inner cell and 
again solving a boundary Riemann problem. Here the problem is complexifled by the fact 
that the limit r — )■ leads to a change of type, from hyperbolic to elliptic, of some of the 
equations. This situation is similar as in the low Mach- number limit, with the additional 
feature that the elliptic equations are one-dimensional, posed along a magnetic field line. 

^For the last reference, see chapters 14 and 15 
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Here, we propose some boundary conditions which are operational and, in particular, 
which avoid the appearance of boundary layers. However, a rigorous mathematical theory 
is not available yet. The geometry of the magnetic field lines plays an important role. We 
assume that the problem is posed in a bounded domain Q with boundary F decomposed 
into 

r = r+ur„uro, (6.13) 

T± = {x eT\ ±b{x) ■ 1^ > 0}, To = {x eT\b{x) ■ 1^ = 0}. (6.14) 

Relative to the magnetic field lines, r_ is the incoming boundary, r+ the outgoing one and 
To, the tangential one. Our proposal of boundary conditions is guided by two principles 
in hierarchical order: first, avoid boundary layers and second, find artificial boundary 
conditions that are as close as possible to the free-space situation (transparent boundary 
conditions). 

The prescription of boundary conditions for n"^ on r_ U r_|_ is guided by the first 
principle, rf satisfies the elliptic equation (16. 9p . The elliptic operator is, at leading order 
in r, equal to — V((6® &)(Vp(n'^) —n^E)). We propose homogeneous Neumann boundary 
conditions associated to the conormal derivative of this elliptic operator on r„ U r_|_, 
namely: 

V\\p{rf) - TfE\\ = 0, on r_ U r+. (6.15) 

In this way, the boundary condition is compatible with the leading order operator inside 
the domain. We can prescribe non-homogeneous Neumann boundary conditions, provided 
that the right-hand side in (16.151) is 0{r). On Fq, we apply the second principle and 
propose a homogeneous Neumann boundary condition: 

u ■ Vrf = 0, on Tq, (6.16) 

which models the flatness of the density proflle near the boundary. This is intended to be 
an approximation of the case where the domain is the entire space. 

For n'^Up we proceed in a similar fashion, considering the elliptic equation (16. 3p . At the 
leading order in r, the principal part of the associated elliptic operator is (6- V)(p'(n'')V ■ 
(fora'^np). The homogeneous Neumann boundary conditions associated to the conormal 
derivative of this elliptic operator on r_ U r_|_, are given by (assuming that p'{n'^) never 
vanishes) : 

V ■ (bn^ul) = 0, on r_ U r+. (6.17) 

The prescription of this boundary condition follows the flrst principle. Again, inhomoge- 
neous boundary conditions can be used provided that the right-hand side of (16.171) is 0(r). 
The second principle leads us to prescribe homogeneous Neumann boundary conditions 
for the parallel momentum on Fq: 

u ■ Vin^ul) = 0, on Tq. (6.18) 
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Now, we consider the transverse momentum and eq. f lG.lOp . The first principle (avoid- 
ance of boundary layers) leads us to propose Dirichlet boundary conditions for n'^u]_: 

ul\r{y) = ulM, Vyer, (6.19) 

with u]_^ satisfying: 

h^uIb = TT^b X {Vp(n^) - n^E} + 0(r), on T. (6.20) 
\Jd\ 

We comment on the conditions (16.151) and (I6.17p . Introducing functions h{n) such 
that h'{n) = p'{n)/n, and 4>{x) such that locally: 

E|| = -V||0, (6.21) 

(I6.15P can be rewritten as 

6 ■ V(/i(n^) + 0) = 0, on r_ U r+. (6.22) 

It expresses that h{n'^) + is locally constant in the direction of the field line on r_ U r+. 
Now, using that V ■ -B = 0, we have 

V-6= -b-\/\n\B\. 

Then, fl6.17p can be recast as 

= 0, on r_ U r+, (6.23) 

which expresses that the quantity n'^uy\B\ is locally constant in the direction of the field 
line on r_ U r+. 

A last comment is that eq. (16. 9p with boundary conditions f l6.15p . (I6.16P or eq. (16. 3p 
with boundary conditions (I6.17p . (I6.18P are well posed for r > but ill-posed in the limit 
r — )■ because the solution is then defined up to a constant (per field line). This induces 
a bad conditioning of these equations when r is small which will require some special 
treatment. We also remark the strong anisotropy of the problems in the direction of the 
field lines. Since h may be time-dependent, we wish to develop solution strategies which 
do not rely on a special set of coordinates related to h. The resolution of these elliptic 
problem will be considered in detail in a forthcoming section. 

7 Time semi-discrete AP scheme 

7.1 Time semi-discrete schemes: general setting 

We rewrite the lEL system using the conservative variables ri^ and = ri^u^: 

dtn^ + V ■ = 0, (7.1) 
rldtq^ + V- r"^/ j I + Vp{n^) = ifE + xB, (7.2) 
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We devise two Asymptotic-Preserving (AP) schemes corresponding to the RIEL-1 and 
RIEL-2 (respectively) reformulations of the Euler-Lorentz system. We start from a dis- 
cretization of system fl7.ip . f l7.2p and design the schemes in such a way that the manipu- 
lations which have led to the RIEL-1 and RIEL-2 reformulations can be performed at the 
discrete level. 

There are several reasons for not using the RIEL-1 or RIEL-2 formulations directly. 
First, these formulations are quite complicated and involve many terms. It is not clear how 
to discretize them in a good way. Second, the scheme must provide consistent solutions 
in both the regimes r = 0(1) and r ^ 1. The RIEL-1 and RIEL-2 forms are adequate 
for the regime r ^ 1 but not for the regime r = 0(1). In this regime, the problem is 
a standard system of conservation laws with source terms, for which a huge literature is 
available (see e.g. jlQl |52l [531 EH]). This literature can be directly adapted to the form 
fl7.ip . (17. 2p but much less obviously to the the RIEL-1 or RIEL-2 forms. For these reasons, 
we develop our schemes starting from (17. ip . (17. 2p . 

We first consider the time semi-discretization, because, in the present example, like 
in many other instances, the design of an AP scheme is primarily a question of time- 
discretization, in a forthcoming section, we will discuss the full discretization of these 
equations by AP methods. The time semi-discretization serves as a preparation for this 
last step. Surprisingly, the algebra is slightly simpler in the discrete than in the continuous 
case. We successively discuss the two reformulations. 



7.2 AP scheme based on the first reformulation 

We start with some notations and preliminaries. We denote by n"^'™, g"^'™ approximations 
of and at time t'" = m5. Since h may depend on time, we index the 'perpendicular' 
and 'parallel' components of a vector field v by the time index m, i.e. 6™ = h{t'^) and 

V = v\\^ir + v^m , v\\^ = vir, v^rn = ir X {v xiT) = iid-w ®ir')v , (7.3) 

and similarly for the parallel gradient and divergence operators. We suppose that E and 
B are known in the course of time and that these projections are available at all times 
without any approximation. The coupling of the Euler equations with a time evolution of 
E and B (e.g. through Maxwell's equations) will not be discussed here, but is of course 
an important and interesting question for future works. 
We now introduce the: 

Definition 7.1 The First Semi-Discrete AP scheme (SDAP-1 scheme) is defined by: 

r {r^(g--^^ - g-) + V ■ f ^""1^" ! ) + 
+V 



n' 

T,m+1 



= ^-."^+1^™+! + g-'^+l X 5™+^ (7.5) 
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The rationale for this approximation is as follows. We start from the following implicit 
scheme: 



= n^."»+i^™+i + X (7.7) 

We then split V ■ into its parallel and transverse components and evaluate the 

transverse component explicitly: 

V ■ q — V||m+i ■ + V_Lm+i ■ 

= V||m+i ■ g^;nt^ + V_l™+i ■ g^C+i + 0{S). (7.8) 

Inserting (17. 8p into (17.61) leads to (I7.4p . Then, we Taylor expand p{n'^'"^'^^) , insert (17.61) 
and use (17. 8p again: 

p(n"''"+^) = p(n"''") + p'(n"''")(n"''"+i - n"'"*) + o(5) 
= p(n"'™)-<5p'(n"'™)V-g"''"+' + o(5) 

= p(n^''") - 6p'{n^n (V|j,n+i ■ gj;::!^ + V^^n+i • q^Zi) + o{6). (7.9) 

Now, we replace p{n'^''^~^^) in (17.71 ) by the right-hand side of (17. 9p . This leads to (17.50 . 

None of the above listed manipulations has altered the conservative character of the 
scheme, nor its consistency. Thanks to these modifications, this scheme is consistent 
with the RIEL-1 formulation of the Euler-Lorentz model, and in the limit r — t- 0, with 
the RIDF-1 formulation of the Drift-Fluid model. This is precisely stated in the two 
following propositions: 

Proposition 7.2 The SDAP-1 scheme can he equivalently formulated as follows: 

ri(n^'™+i - n^'-) + V||™+i ■ gjr+t' + V^^+i ■ glC+i = 0, (7.10) 

T.m+l ^ / T-i ( T.m\ s" // T.m\ (t-i r,m+l i t-i r.m ^ 

^x^+i = X (^V_l™+i ^P[n ' ) - dp (n • ) (V||m+i ■ q^^i^+^ + V^m+i ■ q^^+i^ 



+ r |ri(gl'ri^ - J + (V ■ ( ^ ' ' ))x^.. | J(.7.11) 

rr^gj^rtt' - 5 V||..+i(p'(n-'"^) (Vii^+i ■ gj;::t')) + 5 i^pii vii^+i ■ = 

„r,m „r,m 

= r5 - r (V ■ (^^^^^^^ ))||™+i 

-V||™+i [p(n^''") - 5p'(n^'"^) V^^+i ■ g^-r+i] + 

+ (^-.- _ 5v^™+i ■ glC^O^pi^, (7.12) 

r/ie SDAP-1 scheme is consistent with the RIEL-1 formulation of the Euler-Lorentz model 



54 



Proof: We first take the parallel component of (17.51) and get: 

J 1/ T,m+1 T,m \ , q^'^ ® q^'^ 1 

[q{ra + l -g||™ + l) + (V-( ^^^^^^^ ))||™ + 1^ + 



-Vi 



m + 1 



n^-'^+^Ep+V (7.13) 



Using (17. 4p to eliminate n'^'"*"''^ at the right-hand side of (17.131) . and bringing all terms 
involving q^i^^+i^ to the left-hand side and all other terms to the right-hand side, we find 
fl7l2|) . Now, taking the cross product of ([73]) with 6"^+^ we get flTTT]) . Eq. (Tmj) is a 
time-integrated version of the wave equation (16. 3p . Therefore, the whole SDAP-1 scheme 
is consistent with the RIEL-1 model. This ends the proof. ■ 

Eq. (17.121) takes the form of an elliptic problem for q^^^^^ where the right-hand side 
is known. This elliptic equation, supplemented with the Neumann boundary conditions 
(I6.17p . is well posed and provides the updated value Once is known, rf''^+^ 

can be computed using (I7.10p . Finally, eq. (17. lip , which is clearly a time discretization 
of (16. 2p . can be solved. Alternately (17. lip can be written: 



Id - ^^^^'"^'x j &l = b-^' X (7.14) 

y"'^'^' = (v^-^+i [p{n^n - sp'in^n (Vy^+i ■ gjr+t' + ■ gic+o" 

-n-'-+'E-tl + r |-r^gl'- . + (V ■ ( ^"'"^^/'"^ }) , (7.15) 



where Id denotes the identity matrix and b"^'^^ x denote the matrix of the vector product 
by 6™+^. The vector jg constructed with known quantities and (17.140 can be easily 

solved for g^C+i^ as follows: 



,T,m+l 

m+l • 



This algebraic relation provides the update 
We have the obvious: 

Proposition 7.3 The limit t ^ in the SDAP-1 scheme leads to the following First 
Semi-Discrete Drift- Fluid scheme or SDDF-1 scheme: 

. (7.18) 

-5 V|i™+i(p'(n°'") (Vii^+i ■ gjr+t')) = -Vn^+i [p(n°'™) - 5p'(n°'"^) V^™+i ■ g^C+i] + 

+n°'™+^E"l+\. (7.19) 
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The SDDF-1 scheme is consistent with the RIDF-1 reformulation of the Drift-Fluid model. 



Proof: Eqs. f lTTTD and fmsD are clearly consistent with (ED, (El), while f moll is a 
time-integrated version of (15. 5p . Therefore, the SDDF-1 scheme is consistent with the 
RIDF-1 model. ■ 

This last proposition shows that the SDAP-1 scheme si AP. 

7.3 AP scheme based on the second reformulation 

Definition 7.4 The Second Semi-Discrete AP scheme (SDAP-2 scheme) is defined by: 

ri(n^'™+i - n-'-) + V ■ (gjr+t'6"+' + gl'^+O = 0, (7.20) 
r |5-i(g"'"^+i - g"'") + V ■ ^"'"^f j"'" ^ | + Vp(n"''"+i) = 



The rationale for this scheme is the same as for the SDAP-1 scheme. We start from (17. 6p . 
( 17 ■7p and transform the mass flux using (17. 8p . However, here, we do not transform the 
pressure term in the momentum balance eq. This leads to the SDAP-2 scheme. 

Again, none of these manipulations has altered the conservative character of the 
scheme, nor its consistency. We show that this scheme is consistent with the RIEL-2 
formulation of the Euler-Lorentz model, and in the limit r — )■ 0, with the RIDF-2 formu- 
lation of the Drift-Fluid model. 

Proposition 7.5 The SDAP-2 scheme can be equivalently formulated as follows: 
^§-i^r,m+i _ <5 V||m+i ■ (V||m+ip(n"'™+^) - n"''"+^E™+\) = 

= r6-^n^^- + r |-V ■ q^'"' + 6{V ■ C , (7-22) 

T,m+1 " ^/'^T,ni+lN ^r,m+l rpm+l 

^±•"+1 = 1^^+11 X ( v_L™+ip(ri ) - n ' h^^+j 



(1 ' (x) (1 ' 

+r{5~'W^^ - q2Z.) + (V ■ ))x-+i}), (7.23) 



s-\qtZt' - gfr.O + (V ■ ( ^'J^^ ))r.i + 



+r-i (^V||™+ip(n"''"+^) - n"''"+^Ep+\) = 0. (7.24) 
It is consistent with the RIEL-2 formulation of the Euler-Lorentz model 
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Proof: Taking the parallel component of (17.2 II) . we find f l7.24p . Inserting the value of 
^ll^t^ found from ( \7.2A\\ in (17.201) leads to (17.221) . Taking the cross product of 6"*+^ with 
(17.211) leads to (17.231) . Eq. (17.221) is a time- integrated version of the wave equation (16. 9p . 
Therefore, the whole SDAP-2 scheme is consistent with the RIEL-2 formulation of the 
Euler-Lorentz model. ■ 

Eq. (I7.22P is a nonlinear elliptic equation for the right-hand side of which is known 

from previous time steps. It has a unique solution thanks to the boundary conditions 
(I6.15p . Furthermore, these boundary conditions guarantee that 

V||™+ip(n"'"^+^) - = 0(r). (7.25) 

Eq. (I7.24P looks singular as r — 0. However, with (I7.25p . the seemingly singular term at 
the right-hand side of (I7.24p is of order unity and the equation is not singular as r — >■ 0. 
Finally, (I7.23P can be alternately written: 

(id - 5^^^"^' x) = h^^' X Z-^'/\ (7.26) 



f^T,m ^ T,m 

+r{-r^g-- . + (V ■ 1^ (7.27) 
and has the solution 

^I-i' = ^2^f2|^J+i|2 (-^^'"^'/' + m^'^'l X Z™+V2). (7.28) 

We investigate the limit r — ?■ in the following proposition whose proof is easy and 
left to the reader: 

Proposition 7.6 Taking the limit r — m the SDAP-2 scheme, expanding 72"^'"^+^ = 
^o,m+i _|_ ^j^i,m+i _|_ q(^^ qj^^ using the boundary condition Ii6.15\) leads to the following 
Second Semi-Discrete Drift-Fluid scheme or SDDF-2 scheme: 

V||™+ip(n°''"+') - n°'™+^E™+^i = 0, (7.29) 
^-i^o,m+i _ 5 v,|™+i ■ (V||™+i(p'(n°'"^+i)ni''"+i) - n^'^+iEp+iJ = 

= 5-'n''- + -V ■ + 5(V ■ ( ^ „J ))r.i , (7.30) 



„0,m ^ „0,m, 

c-_l/ 0,m+l 0,m \ , /-r-r /'i ^ 'i w , 

d -g|i™+i) + (V- ( — — ))ii-+i + 



(^V||-+i(p'(n°'"^+^)ni''^+i) -n^'^+iElp+l) =0. (7.31) 



q^t' = X (V^™.ip(n°'-+^) - n°'™+iE7+,\), (7.32) 
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The SDDF-2 scheme is consistent with the RIDF-2 reformulation of the Drift-Fluid model. 

This last proposition shows that the SDAP-2 scheme si AP. 
7.4 Time semi-discrete schemes: conclusions 

In the previous sections, we have derived two different semi-discrete AP schemes (the 
SDAP-1 and SDAP-2 schemes). In the limit r — ?■ 0, these two schemes are respectively 
consistent with the two previously established reformulations of the Euler-Lorentz problem 
(the RIEL-1 and RIEL-2 reformulations). To derive these schemes, we have started from 
the time continuous problem in its original form (the lEL form) instead of using the 
reformulated forms. Both schemes are derived from an implicit scheme where the flux in 
the mass conservation equation, the pressure flux in the momentum conservation equation, 
and the Lorentz force are evaluated implicitly. Then, the two schemes are transformed 
somehow similarly: in the SDAP-1 scheme, the implicit pressure gives rise to an elliptic 
equation for the parallel momentum through the use of the mass conservation equation. 
A symmetric operation is performed on the SDAP-2 scheme, where the implicit mass 
flux gives rise to an elliptic equation on the pressure through the use of the momentum 
balance equation. Both elliptic equations are degenerate: they are one-dimensional elliptic 
equations posed in the direction of the magnetic field lines. 

The derivation of time semi-discrete schemes is a preparation for the development of 
the fully discrete ones which will be performed in the next section. Time semi-discrete 
schemes are also convenient for linearized stability analyses in the spirit of section 11.41 or 
12.21 The stability analysis of these schemes will be considered in future work. 

8 Fully discrete AP scheme 

8.1 Fully discrete schemes: general setting 
8.1.1 Classical schemes 

Before considering AP schemes, we recall the framework of the classical explicit shock- 
capturing schemes. We start from the original lEL formulation of the problem ( 17. ip . (17. 2p . 
which we write as follows: 



3 




(8.1) 
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where U = U{x, t) is the vector of conservative variables, fj{U), the flux and g, the source 
term: 

U=( fj{U) = f _i , , V (8.2) 

We denote by Cj the unit vector in the j-th direction (for instance, 62 = (0, 1,0)). Here, 
we always assume that the space dimension is 3. 

We develop a finite-volume formulation on a structured, cartesian mesh. However, the 
concepts would easily be generalized to a finite volume method on an unstructured mesh. 
Let K = {ki,k2,kz) G be a multi-index, and Ck = hi[ki — 1/2, ki + 1/2] x h2[k2 — 
1/2, k2 + 1/2] X h-^lk^ — 1/2, k^ + 1/2] be the associated finite- volume cell, with space steps 
hj in the j-th direction. We denote by f/|^ an approximation of U{xK,t"^), where xk = 
{kihi, k2h2, k^h^) is the center of cell ck- Similarly, /j |^_,_e^/2 denotes an approximation 
of fj{U{xK+ej/2, ^'")), with XK+ei/2 = ((^1 + l/2)/ii, A;2/i2, ^s^s), and similarly for j = 2, 3. 
Finally, g'^iU) ^ g{t"',U). 

The classical schemes are written as follows 

Definition 8.1 Classical explicit shock- capturing explicit schemes are given by: 

3 

The fluxes are the sum of a central discretization term and a numerical viscosity term, 
according to: 

fj\K+e,/2 = \ [fmK) + fmK+e,) + ^^^K+e, /2{U\k ' f^l^+e,)] , (8-5) 

with a suitable viscosity matrix /^j|^+e ./2- 
We denote by 

the numerical viscosity. In explicit shock capturing methods, /^j|^+e./2 is derived from 
the Jacobian matrix of the flux functions. For instance, the Roe scheme corredponds to 



M'j\K+ej/2 



Of J 

-gfjiU\K+e,/2] 



•7) 



where f^|^+ej72 ^ conveniently chosen average state between f/|^ and U\]^^^_.. The 
Rusanov scheme would correspond to /^j|^_|_e^/2 begin a scalar such that 

H\K+e,/2 = max{max{ \XMk)\ , |A,(f/|^+e,/2)l , |A,(t/|x+e,)l }, 

df- 

such that Xj{U) is an eigenvalue of -^{U)}, (8.8) 
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The CFL stability condition, which guarantees the stabihty of the scheme, is as follows 



ji5 <h = min/ij, (8-9) 

with 

11= max 1x^+6/2 5 (8.10) 

Any other shock capturing methods can be considered as well. We refer the reader to 
PUI [521 l53l [63] and references therein. 

8.1.2 AP schemes 

By contrast to our presentation of the Semi-Discrete schemes in section [TJ we gradually 
derive the final expression of the fully-Discrete AP schemes from a general semi-implicit 
shock capturing scheme framework. In this section, we present the common starting 
point for the two AP schemes which we consider in this work. Then, in two forthcoming 
sections, we will develop the specificities of each of these schemes which are intended to be 
consistent discretizations of the Reformulated Euler-Lorentz models RIEL-1 and RIEL-2. 

We plan to design AP-schemes from minor modifications from classical shock capturing 
schemes in order to ensure that the discretization of the left-hand side of (17. ip . (17. 2p is 
in conservative form. The conservativity property guarantees correct shock speeds at the 
discrete level. It cannot not be guaranteed if the scheme is developed from the RIEL-1 
or RIEL-2 reformulated systems because of the presence of many and sometimes complex 
terms in these formulations. Another reason for dealing with the original system (17. ip . 
(17. 2p is the need for numerical viscosity to stabilize the discretization. While it is easy 
to adapt the literature to (17.11) . (17. 2p . it is uneasy to decide where and how numerical 
viscosity should be added to the RIEL-1 or RIEL-2 formulations. 

The common framework for both AP schemes is a modification of (18.41) where some 
kind of implicitness in the flux and source terms is introduced: 

3 

6~\U\T' - U\-) + Y.hj\WXl/2 - UT-eJ = r-'g'-^'m-^^ (8.11) 

where the implicit fluxes are denoted /j |^+ey2 distinguish them from the explicit ones 
(18. 5p . The AP schemes will be such that the implicit fluxes /j|^+e./2 fairly simple 
modifications of the explicit fluxes (18. 5p . Like in the time semi-discrete case, we construct 
them by making the mass and pressure fluxes implicit. Additionally, the implicitness only 
applies to the central part of the flux, because impliciting the viscosity is not needed to 
make the scheme AP. The implicit fluxes are thus given by: 



im+l 



1 



fj\K+ej/2 2 



f^T'+f^TA+^^^K^eAu\K-u\'^^^) . (8.12) 
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The bars denote implicit central fluxes fj\x given by: 

r \m+l 



JJ\K 



im+l 
'ijlK 



(n-ig,g)|^ + r-yn|r')e, 



(8.13) 



and {n~^qjq)\']^ is a short-hand writing for (^~^)|^ qj\^ q\^- Note that this part of the 
momentum flux is kept explicit, because impliciting it is not needed to make the scheme 
AP. The viscosity matrix and the CFL conditions are constructed only from the explicit 
part of the system, i.e. they are associated to the flux: 





n~^qjq 



(8.14) 



This systems has eigenvalues and uj = qj/n. For instance, for the Rusanov scheme, 
A^il^+e /2 i^ ^ scalar equal to the maximal value of \u\ in the adjacent cells: 



A''jlx+e,/2 — niax{ \u\^\ , 



ej/2\ 5 Wk+s 



}• 



(8.15) 



In doing so, neither the numerical viscosity nor the CFL condition depend on r, a condition 
for the scheme to be AP. Any other shock capturing methods can be considered as well. It 
may improve the stability of the scheme to keep some small explicit part in the mass and 
pressure fluxes, in the spirit of the method proposed in [30]. We will defer the development 
of this idea to future work. 

Inserting (18. 6 p and (I8.13P into (I8.12p and the resulting expression into (18. 4p . we get 
the following expression for the scheme: 



1/ im+l 



[n 



K 



m+l 
n]\K+ej/2 



f I m+l \ 

Jn'i\K-ej/2) 



0, 



(8.16) 



3 



I X ) + ^ ^ ( fqd I X+ej /2 fn 



m+l \ 
<li]\K~ej/2> 



— 1/ im+l 7-1 im+l 



r [n 



K 



|m+l^ 



1,2,3, 



^.17) 



Here, we have supposed that the electric and magnetic fields are appxomitated by cell- 
centered discretizations E\^, B\^. The fluxes are given by: 



m+l _ _ 



I m+l I I m+l I r-) |m 
^j\k ~^ 1j\ K+ej ^nj \ K+Cj /2 



I m+l 



r I m+l 

J<}d\K+ej/2 ~ 2 

1 

+ - 



|m+l^ 

'j\k ) 



PKli^+e, 



im+l 



5jl + 



^,j = 1,2,3, (8.19) 



where D^j l^+e 12 ^^"^ ^Qijlx+e 12 entries of the numerical viscosity vector -Dj |^ 



3\K+ej/2 ■- 



Dr. 



+e,/2- 



(8.20) 
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and D 



nK+ej/2 



is given by fl8.6p (but with the viscosity matrix associated to the exphcit 



part of the flux, as pointed out in the previous paragraph). For later usage, we define: 
1 



^ qij\K+ej/2 — 2 



z,j = 1,2,3. (8.21) 



Alternately, inserting iKW^ . iKW^ into f lgl^ . (ISTr]) . we can write: 

2 



„ |m+l im+l 



_l_ A I'Ti _ n 



m 

qi\K 



r-^(nir'i?.ir' + {q\T' X ^ = 1,2,3, 



with 



2 

3 



A I™ 



A 



n .1™ _ n .I'" 



jp \m TP \m 

^qij\K+ej/2 ^qij\K-ej/2 



i = 1,2,3. 



(8.22) 
(8.23) 

(8.24) 
(8.25) 



We are now going to modify this scheme in two different ways in order to make each of 
the resulting scheme consistent with either the RlEL-1 or the RIEL-2 formulation of the 
Euler-Lorentz model, and a fully discrete counterpart of the time semi-discrete SDAP-1 
and SDAP-2 schemes respectively. 



8.2 Fully discrete AP scheme based on the first reformulation 

First, we decompose the momentum into its parallel and perpendicular parts. Specifically, 
we denote by 



(Id - 6 ® 6)1™ = Id - 6|™ ® g^l^ = (Id - 6 ® 6)|^g 

qwiK - qiK ■ o\k, qiK - q\\\KO\K + quK- 



(8.26) 
(8.27) 



Therefore, we can write 



im+l 

q\K 



q\\\T^b\T' + (Id-6®6)|^+'g|^+\ 
q\\\T'b\T' + (Id-6®6)irMx + 0(5) 



(8.28) 



and insert this approximation into (18.181) . This leads to a modified expression of the 
discrete mass flux: 



f im+l _ ^ 
Jnj\K+ej/2 ~ ^ 



|m+l7 im+l I im+l 



+ D 



nj\K+ej/2^ 



(8.29) 
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with 



D. 



nj\K+ej/2 — ^ 



+ 



(8.30) 



The discrete mass balance eq. f l8.16p with the modified fiux fl8.29p is now written: 



K 



— n\ 



K. 



m+1 



+ An\T< = 0, (8.31) 



with 



A 



n\K 



D 



nj\K+ej/2 ^nj\K-ej/2 



(8.32) 



Proceeding hke in section [7]2l We now expand p{nj\'^'^^), using f l8.22p and fl8.3ip : 



+ 



i'i^-5y(ni^)X^^U 1^:^^,6,1: 



■m+l 



(8.33) 



+A„|^ + 0(5)), (8.34) 
<l\^T\MT'e'\^0{b\ (8.35) 



with 



^l^=P(^lx)-^PVlA-)An|^. 



(8.36) 



Then, replacing p(n|^"''^) in fl8.19p by its approximation f l8.35p leads to a modified mo- 
mentum flux: 



Im+l p |m 



+ 



r=l 



K+e 



E r I m+l 

9 mK+ 



■1' ^ 2 

r=l 



I m+l 



m+l 



|m+l 



e,+er "''lA'+e,+er ^11 lif+e,-er "^l/f+e,-e. 



5.,, (8.37) 



with 



F 



I m 

lE'+ej/2 



r 



-P|^ + -P|^+e,' + -^'?iilx+e,/2) "^^i — 1) 2, 3. 



(8.38) 
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The momentum balance eq. (18 .17^ with the modified fiux fl8.37p is now written: 

/)~^('/7-l'"+^ — -I- A I™ 



U-1 



„ |m+l L \m+l _«,,r'"+l h \m+l 



i=i 



im+l L im+i im+i l 

y|| lii:-ei+ej "ilX-e.+e,- yill/f-ej-ew " 



|m+l 



m+1 



|m+l 



r-i(n|-+^E,|™+i + x B\^+%), i = 1,2,3, 



with 



3 



P im p I'm 

^ qij\K+ej/2 ^ qij\K-ej/2 



, z = l,2,3. 



^.39) 



(8.40) 



Thus, we are led to the following definition: 
Definition 8.2 The first Fully-Discrete AP scheme (or FDAP-1 scheme) is given by 



1/ im+l 



K 



n\ 



7 I m+1 _„,,|rrt+l t im+l 
y||lii:+e, '^i\K+ei '^\\\K-e, "jlK-Ci 



A„|™ = 0, (8.41) 



A I™ 



„ im+l L |m+l 7) 



'd\\\K-ei+ei "jlK-e^+Cj 'd\\\K-e,-ej 



m+1 



I m+1 



E"7' 

= r-i(n|-+iE,|-+i + (g|-+i x 51^+^),), z = 1,2,3, (8.42) 



wt/i A„|^ (/zwen &y / fO^) . / fO^) . / fO) on i/ie one /ianfi and A^J^ {^^^j, / f05l) . 
/ fO^) . ( fOTl) . / fO) on i/ie oi/ier /iond. 



The quantities A„|^ and AgJ^ only depend on known quantities at the previous time 
steps (with the exception of the magnetic field which is supposed to be know but taken 
at the current time step). A„|^ collects a consistent approximation of V_l ■ q± and 
the contribution of the numerical viscosity (see (18.301) ). The quantity Aq.|^ collects 
a consistent approximation of V ■ {n~^q ^ q) + r^^Vp(n) and the contribution of the 
numerical viscosity (see (18.381) and (18.211) ). 

We now show that this scheme can be recast into a consistent approximation of the 
RIEL-1 reformulation and that it is actually AP. 
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1/ im+l 



6-\n 



K 



A 



Proposition 8.3 The FDAP-1 scheme can he equivalently formulated: 
2-^ A 



771+1 

K 



I m+l 



\m+l 



n\K 



(8.43) 



m+l 



I m+l 



K+Si—Cj j\K+ei~ej 



«J=1 



q\\\K-e,+ej " 



I m+l 



„ I m+l L 

Ei+ej "Jl_ft:-ei+ej 'i\\\K-ei-ej "j 



I r2 7-1 I m+l \ ^ J ^ im+l t im+l ^ |m+l t |m- 



im+l 



i=i 



q± 



m+l 



5 LB| 



m+l 



m+l I "IK 
K 



X gj 



I m+l 



^1 



m+l 
K 



B 



m+l I 
K 



{G\]^-n\T'E\r'-r5~'q\]^}. 



with 



(8.44) 



^.45) 



i?|- = i?|- + 5E||ir^(n|--5A„|™), E| 

3 3 



I m+l 
IK 



1=1 



R\K = rYl b^\T' Qi\K -rSYl b^\T'^^ 



m 
qi\K- 



i=l 

—lf~i |m X |m 

^ <-^i|_fi' — ^(ji \k 



i=l 



U-1 



ejX^^i^ ?||IS^+ei+ej^ ^ilS^+e,+ej 9|| 



(8.46) 
(8.47) 



m+l 1 I m+l 
\K+ei-ej ^jlK+ei-Cj 



i=i 



I m+l 



I m+l 



I m+l 



(8.48) 



i=i 



This scheme is consistent with the RIEL-1 formulation of the Euler-Lorentz model and is 
a full discretization of the SDAP-1 time semi-discrete scheme. 

Proof: We first consider tlie parallel component of the momentum and take the dot 
product of (18.421) by (i.e. we multiply (I8.42p by ftj|^+^ and sum over i). We 

multiply by t5 and we find: 

1^ 4 

«,i=i 



m+l 
K 



|m+l 



A'+ei+e,' iK+ei+e,^ 



511 1 



m+l 



I m+l 



iC+ei— e,- jl/C+ei-e; 



-nVnl™- \h |™+1 n h n A |"l+l 1 

y \"'\K-ei) "Ak \y\\\K-ei+ej '^jlK-ei+Cj 'd\\\ K -ei-ej '^j \ K-ei-ej j 



Sj2^\T'h\T'Ei\T' + R\k^ (8.49) 



i=l 
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with given by f l8.47p . Inserting f l8.4ip into (18 .49^ allows us to replace n\]^~^'^ by an 
expression involving and known values at time n. This leads to f l8.44p . We now 

consider the transverse component of the momentum. Taking the vector product of f l8.42p 
with and dividing by we find (gM>- 

We now show that the FDAP-1 scheme is consistent with the SDAP-1 time semi- 
discrete scheme (and consequently, with the RIEL-1 reformulation of the Euler-Lorentz 
model). Eq. (I8.44p is a discrete version of the elliptic equation: 

rgy -6H- {V{p{n)V ■ (gyfe)) - EV ■ (g||6)) = R. (8.50) 
Developing the term R\]^ using fl8:iOD . dHjH]), fl836D . dHJH), we find: 

c 3 



2=1 



n\K) 



6 



i=l j=l 

3 3 



2E^^IS^' E^^' [D,a\K^e,/2-D,M-e,/2\ + Y.^^\T' Q^\K \ ■ (8-51) 
i=l j=l i=l J 

We easily convince ourselves that 6~^R\]^ is a consistent approximation of the right-hand 
side of (I7.12p . and consequently, that (I8.44p is consistent with (I7.12p . Eq. (I8.43P is clearly 
a consistent approximation of f IT.lUp . Finally, appears as a consistent approximation 
of Vp(n) -|- rV ■ {n~^q ® q) and therefore, fl8.45p is a consistent approximation of (17. lip . 
This ends the proof. ■ 

Supplemented with the Neumann boundary conditions (16.170 . the elliptic equation (I8.44p 
is well-posed. Its inversion allows us to compute the values of gy l^"*"^ from the right-hand 
side R\]^ which only involves known values. From the knowledge of ^'nl^'''^, we can use 
(I8.43P to find the new values n\^'^^ of the density. Finally, the quantities lying at the 
right-hand side of ( 18.45P are known at this level of the recursion. Therefore, can 
be computed thanks to the formula: 



6 B 



m+l I 



T 



2 



+6 6|^+'x) {G|^-n|^+^E|7^+i-rr^g|^}, (8.52) 



where the right-hand side is known. 
We now investigate the limit r — )■ 0. 
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Proposition 8.4 The limit r — o/ the FDAP-1 scheme is the following First Fully- 
Discrete Drift-Fluid scheme (FDDF-1 scheme): 



r'(nir' - "IS) + E 



2 



I m+1 L I m+1 I m+1 l |m+l 

y||lii:+e,- "jlX+e,- ^\K-e, '^j\K-e. 



A 1™ — n 



.53) 



y|ilA'+ei+ej "jl_ft:+ei+ej Vll I _ft:+ei-ej "jlx+ 



m.+l 



im+l 



|..-,^ L im+l im+l 7 im+l 

'i\\\K-ei+ej "j\K-ei+ej H\\\K-ei-ej "j\K-ei~ej 



+ 



+^ ^\\\k ^l|li^+ej ''ili^+e, - ^ 



I m.+l L I m+1 



„ im+l 

Q±\k 



m+1 / I m 



n 



K 



5An 



i=l 



m+1 
K 



:^x{GX-n\T'E\^+'}. 



.54) 
.55) 



with 



Go I m 
i \k 



m 

K+e, 



6h7^ 



p\m 1 



I m+1 
y||lit+e,+ei 



I m+1 



im+l 



ei+Gj "J lit+Ci+ej y|| li^+e—ej lit+Ci-ej 



I m+1 



m+1 



I m+1 



m+1 



.56) 



This scheme is consistent with the RIDF-1 formulation of the Drift-Fluid model and is a 
consistent space- discretization of the SDDF-1 time semi-discrete scheme. 



The last two propositions show that the FDAP-1 scheme is AP. We now consider a scheme 
based on the second reformulation. 



8.3 Fully discrete AP scheme based on the second reformulation 

Following the strategy developed in section I7.2[ the second Fully-Discrete AP method 
(or FDAP-2 scheme) consists in using the modified mass fiux fl8.29p but the unmodified 
momentum fiux fl8.19|] . Therefore, the FDAP-2 scheme is defined by: 
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Definition 8.5 The Second Fully-Discrete AP scheme (or FDAP-2 scheme) is given by 



6'\n 



1( ^ |m+l 
K 



n\ 



im+l L im+l 



m+1 7 im+l 



^ h-^ r 



A„|^ = 0, (8.57) 



1,2,3, 



58) 



with A„|^ given by i\8. Ji^j) . ([K^j o"- ^/^e one /ianc? and A^. |^ zs g'zt'en 6?/ /[8.25\) . 

i\8.21\) on the other hand and we recall that the definition of the numerical viscosities is 
given by l{8.6\) . 

The quantities A„|^ and Ag, |^ only depend on known quantities at the previous time 
steps. A„|^ collects a consistent approximation of ■ q± and the contribution of the 
numerical viscosity. The quantity A^. |^ collects a consistent approximation of V-(n~^g® 
q) and the contribution of the numerical viscosity. We now show that this scheme can be 
recast into a consistent approximation of the RIEL-1 reformulation and that it is actually 
AP. 

Proposition 8.6 The FDAP-2 scheme can be equivalently formulated: 



^ 2 



' im+l I m+1 



m+1 ^ 

\K+ej+ei/ 



— p{n 



m+1 



— n\ 



m+1 rp |m+l 
X+e,-^«lK+e, 



-b 



m+1 7 I m+1 
'j\K-ei '^i\K~e, 
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m.+l ^ 
K—ej+Ci ) 



n 



+r5-V|2+^ + rA 



m+1 77' I m+1 
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i=i ^ 



n\K 
m+1 
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6 B 



^|m+l 



im+l 

\k 



\B 



m+1 1 
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X \G 



\K 



n 



m+1 j^\m+l 



K 
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with 



^ h-^ r 



m+1 



im ^ 



|m+l 



Q\K~e,'^j\K-ei 



i=l 



-1^ im 



U-1 



W\T+\:) - p{n\T.\)] + K\k- 
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(8.59) 
(8.60) 
(8.61) 

(8.62) 

(8.63) 
(8.64) 



This scheme is consistent with the RIEL-2 formulation of the Euler-Lorentz model and is 
a full discretization of the SDAP-2 time semi-discrete scheme. 



Proof: Taking the parallel component of (18.581) . we get f l8.60p . Inserting it into f l8.57p 
leads to (18 .59 p . We now consider the transverse momentum. Comparing (I8.58P with 
(I8.42P and having the expression (I8.48P of in mind, we notice that relation (18. 45 p . 
which was derived in the case of the FDAP-1 scheme, applies to the FDAP-2 scheme as 
well with simply replaced by This leads to (I8.6ip . 

We now show that the FDAP-2 scheme is consistent with the SDAP-2 time semi- 
discrete scheme (and consequently, with the RIEL-2 reformulation of the Euler-Lorentz 
model). Eq. (I8.59P is a consistent approximation of 

-5V ■ {{h ® h){Vp - uE)) + tS-^u + tA = 0. 

Furthermore, a more careful inspection of A„|^ easily shows that it is a consistent ap- 
proximation of the right-hand side of (17.221) . Therefore, Eq. (I8.59P itself is consistent with 
(17:22|) . It is then clear that (IHIeOj) is consistent with ([7211). Then, the quantity is 
a consistent approximation of Vp(n) + rV ■ {n^^q q) and therefore, (18.610 is consistent 
with (17:231) . I 

Eq. (18.59P is a discrete elliptic equation. Supplemented with the Neumann boundary 
conditions (I6.15p . it is well-posed. Its inversion allows us to compute n\^'^^ from the 

quantity A„|^ which only involves known values. Once n\^'^^ has been found by the 
inversion of (18.591) . the parallel momentum eq. (I8.60p can be used to find Q'ljl^'''^. Finally, 
the quantity is known from the previous steps of the recursion. Eq. (I8.52p applies 
with Gi\^ replaced by and determines qiXx^^ from the known values appearing at 
its right-hand side. 

We now investigate the limit r — 0. 

Proposition 8.7 Taking the limit r ^ in the FDAP-2 scheme, expanding n'^l'^'^'^ = 
^0|m+i _|_ ^^i|m+i _|_ ^^^^ (where the dependence of the solution upon r has been restored) 
and using the boundary condition Ii6.15\) leads to the following Second Fully- Discrete Drift- 
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Fluid scheme or FDDF-2 scheme: 
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r/iis scheme is consistent with the RIDF-2 formulation of the Drift-Fluid model and is a 
consistent space- discretization of the SDDF-2 time semi-discrete scheme. 



Proof: Expanding (I8.59P in powers of r, we find at the leading order: 
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and at the next order, fl8.66p . With a discretized version of the homogeneous Neumann 
boundary conditions fl6.15p . eq. (I8.69P can be integrated once and leads to fl8.65p . Taking 
the r — )■ limit in fl8.60p leads to f l8.67p . because the leading order (in r~^) term cancels 
due to (18.651) . Finally taking the r — )■ limit in fl8.6ip leads to f l8.68p . This scheme is obvi- 
ously consistent with the SDAP-2 scheme and consequently, with the RIDF-2 formulation 
of the Drift-Fluid model. ■ 



The last two propositions show that the FDAP-2 scheme is AP. 
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8.4 Fully discrete AP scheme for the Euler-Lorentz model: con- 
clusion 

In this section, we have derived two fully-discrete schemes for the Euler-Lorentz model, 
which are consistent with the Drift-Fluid limit when the parameter r (representing the 
scaled cyclotron period and Mach number) tends to zero. They have been shown to 
be respectively discretizations of the two reformulations of the Euler-Lorentz model, the 
RIDf-1 and RIFD-2 models. These schemes are based on standard explicit shock-capturing 
schemes where implicit evaluations of the mass and pressure fluxes and of the source terms 
have been introduced. Then, some simple approximations have been performed. These 
approximations alter neither the conservative character of the scheme, nor its consistency 
but give rise to discrete elliptic equations for the parallel momentum and the density 
respectively. These elliptic equations are strongly anisotropic, with a diffusion which is 
concentrated along the magnetic field linres. Therefore, they lead to AP schemes provided 
that they are uniformly well-posed as the parameter r tends to 0. The uniform invertibility 
of these elliptic equations as r tends to as well as a practical methodology to solve them 
in coordinate systems which are independent of the magnetic field lines are investigated 
in the next section. 

9 Numerical resolution of strongly anisotropic ellip- 
tic problems 

9.1 Introduction to strongly anisotropic elliptic problems 

This section is concerned with the numerical invesion of strongly anisotrpic discrete elliptic 
problems such as those undelying the FDAP-1 scheme (18.44^ or the FDAP-2 scheme (18.59p . 
More specifically, we are interested in: 

Definition 9.1 The continuous anisotropic elliptic (AE) problem is as follows: findn'^{x) 
defined on Q as the solution of 

Tn7 - V|| ■ (Vyn^ - E\\) = tF, m fi, (9.1) 

with boundary conditions: 

{h-v){V\\n^ -TfE\\) = Tg, on r+ur_. (9.2) 

Here, the domain geometry is as defined in section \6.3[ The vector field E{x) and the 
function F are known smooth functions defined on fl while g is a known function defined 
on r_|_ U r_ . The parameter r is positive. 

This problem is the generic one which needs to be solved when using the SDAP-2 scheme, 
specifically when inverting fl7.22p . with boundary conditions (16.151) up to obvious nota- 
tional changes. Indeed, replacing by r at the left-hand side of (17.221) . and lumping all 
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the terms at the right-hand side of f l7.22p into F, we find (19. ip . We also have to assume an 
isothermal pressure relation p{n) = n in order to get the linear AE problem. A nonlinear 
pressure relation would give rise to a nonlinear anisotropic elliptic problem: 



The resolution of this nonlinear problem by means of Newton's method would lead to 
solving a sequence of linear problems of the AE type. It may be desirable to design more 
elaborate methods in the small r case but these developments will be the subject of future 
work. In these notes, we will restrict ourselves to the linear AE problem. 

The AE problem enters the class of strongly anisotropic problems because it can be 
recast in the form 



Thus, it appears as an elliptic problem with a diffusion matrix equal to b b. This 
matrix is a rank one matrix, and is therefore singular on a two-dimensional manifold (the 
orthogonal plane to b). This is highly degenerate elliptic problem. 

The AE problem is also nothing but a one-dimensional elliptic problem posed along the 
field lines. However, we aim at designing methods which do not require the computation 
of these field lines, nor any integration along them. The difficulty with the AE problem is 
that it becomes singular when r — )■ 0. Indeed, letting r — )■ in the AE problem leads to 
V||n° — n^E\\ = 0. Supposing E = just for the sake of simplicity, this condition becomes 
V||n° = and means that n° is constant along the magnetic field lines. However, this 
constant is undetermined by the leading order equations in the AE problem. To find the 
value of this constant, it is necessary to expand n'^ in powers of r: n'^ = + rn^ = o(r) 
and to look for the existence condition for the first order perturbation n^. From the 
numerical viewpoint, any standard discretization of the AE problem (like finite-difference 
or finite-element methods) will lead to inverting a matrix with condition number of the 
order 0(r~^), and is therefore unpractical if r ^ 1. The goal of the following discussion is 
to propose inversion methods with uniformly bounded condition number with respect to 
r when r — )■ 0. We will again call these methods 'Asymptotic-Preserving' (AP) methods 
for highly anisotropic elliptic problems. 

These methods can be extended when the elliptic operator has a transverse part of 
order 0(1) (when r — )■ 0). Such problems are encountered for instance in ionosphere 
models (see []) and the presented material has been the subject of a []. The remainder 
of this section is organized as follows: we will first discuss a simple one- dimensional 
framework in the continuous case. Then, we extend it to the general three-dimensional 
framework. The, we discuss the discrete case, specifically aiming at solving the set of 
equations resulting to the application of the FDAP-2 scheme for the Euler-Lorentz model. 
Again, we start with a toy one- dimensional problem and then, in the last section, we 
discuss the fully discrete Three dimensional problem arising from the FDAP-2 scheme. 



rn^ - V|| ■ (y\\p{n^) - n^E\\) = tE, in Q, 
{b ■ z/)(V||p(n^) - n^^ii) = rg, on T+ U r_. 



(9.3) 
(9.4) 



rn^ -V ■ {{b®b){yrf -rfE)) = tE, in 
{b ■ u)b ■ {Vrf - n^E) = rg, on r+ U r_. 



(9.5) 
(9.6) 
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9.2 A simple one-dimensional "anisotropic" elliptic problem: 
the continuous case 



This section is intended for introductory and illustration purposes. Basically, this one- 
dimensional example is illustrating what happens along each magnetic field line. The 
ID-AE problem is defined as follows: 

Definition 9.2 The one- dimensional continuous anisotropic elliptic (ID-AE) problem is 
as follows: find n'^{x) defined on [0, 1] as the solution of 

Tn" - i^rf - rfE] = tF, m [0, 1], (9.7) 
dx \dx J 

with boundary conditions: 

-^n7 — ri^E = rg, at x = and x = 1. (9.8) 
dx 

The vector field E{x) and the function F are known smooth functions defined on [0, 1] 
while g is a known function defined only at x = and x = 1. The parameter t is positive. 



The first order derivative term can be removed by means of a simple transformation. 
Introducing the electric potential 0(x) such that 

(j){x) = - I E{y)dy, (9.9) 
Jo 

we remark that 

d diT 

-n^E = e-^-^, = e^n\ (9.10) 
dx dx 

Then, the IDAE problem is written as follows: 

Proposition 9.3 With defined by Ii9.10\) . problem IDAE is equivalent to the following 
one-dimensional modified anisotropic elliptic (ID-MAE) problem: 

re-^u^ - 4- (e"^^] = rF, m [0, 1], (9.11) 



dx \ dx 

idu^ 



dx 



rg, at X = and a; = 1. (9.12) 



For the limit r — )■ of the ID-MAE problem, we have: 
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Proposition 9.4 Let he the solution of the ID-MAE problem. Then, when r — )■ 0, 
converges to where is constant on [0, 1] and given by: 

• 1 

F{x) dx + [g]l 

= ^1 , (9-13) 

e-'^(^) dy 



where [g]l = g{l) - g{0). 

Proof: We expand u"^ = n° + ru^ + o(r). Inserting this expansion in the ID-MAE 
problem, we get, at leading order: 

4:(e-'^)=0^ in [0,1], (9.14) 



di'tij V di^jC 
du^ 

e-^— =0, at a; = and X = 1. (9.15) 

dx 

and at the next order: 

e"^— = g, at X = and X = 1. (9.17) 
dx 

From fl9.14p . fl9.15p . we get that is a constant over [0, 1]. Then, for problem fl9.16p . 
(19.17P to have a solution, a compatibility condition is required. Indeed, integrating ( 19.16^ 
upon X G [0, 1] and using the boundary conditions (19.171) leads to condition (19.131) . ■ 

Back to the n'^ variable, we find that n'^ — )■ ra" with = vPe~'f' with the constant vP given 
by (I9I3]). 

We now introduce a variational formulation. We denote by V the space i^^(0, 1) = 
{u e L^(0, 1) \ u' G L^(0, 1)} where L^(0, 1) is the space of square integrable functions 
on [0, 1] and u' denotes the space derivative of u. We endow V with the norm \u\v 



\u\ + \u 



/p)i/2 where |n| = {^^ \u[x)\'^dxY^'^ is the norm on L^(0, 1). The scalar products 
on V and L^(0, 1) are respectively denoted by {u,v)y and {u,v). We also define the 
following bilinear forms on V and //^(0, 1) respectively: 

a^{u,v) = I e~^^^dx, {u,v)^= f e~'^uvdx. (9.18) 
Jo dx dx Jq 

Finally, we define the linear form F on ^ by: 

{T,v) = {F,v) + [gv]l (9.19) 

Now, taking v E V , multiplying f l9.1ip by f , integrating on a; G [0, 1] and using the 
boundary condition (I9.12p . we are led to the following 
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Proposition 9.5 A variational formulation of problem ID-MAE is as follows: 

Find such that: 

a^{u\v) + T{u\v)^ = T{V,v), "iveV. (9.20) 



We note that the bihnear form a^{u, v) + r(u, v)^ is coercive on V . However, the best 
coercivity estimate is of the order of r: 

a^{u,u) + T{u,u)^>CT\v\y. (9.21) 

This leads to an estimate of the solution in terms of the right-hand side F as follows 

C 

W\v<-\^\v'- (9.22) 
r 

This estimate deteriorates as r — t- 0. Therefore, the condition number of any standard 
finite element method based on (I9.20p tends to oo as r — 0. Such a method cannot 
be used when r <C 1. We aim at finding a variational formulation which has bounded 
condition number when r — > 0. For this purpose, we adopt the method of [IT] based on 
a micro- macro decomposition of the solution . We indeed view the limit of when 
r — 7- as the macroscopic component of and the difference — vP as a, microscopic 
correction. To define the functional spaces associated to the macro and micro components 
live, we introduce the following decomposition: 

Definition 9.6 We define the spaces G and A as follows: 

G = {v eV\v' = Q}, A = {v eV\v{Q) = Q}. (9.23) 

Proposition 9.7 (i) We have: 

V = G®A. (9.24) 

(a) G is characterized by 

ueG ueV such that a^{u,v)=0, \/v e A. (9.25) 



Proof: (i) is obvious, 
(ii): that u E G =^ a^{u,v) = 0, Wv G 
Conversely, suppose that u E V such that 
with p E G and q E A. Then, = a^{u,v) 



A is obvious (it is even true for all v E V). 

a^{u,v) = 0, Wv E A. Decompose u = p + q, 
= a^{q, v), Wv E A. Choosing v = q, we have 
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(^4>il^^) — 0- -B^t' by Poincare inequality, y/a^ is a norm on v and is equivalent to | ■ |y. 
It follows that g = and that u = p E G. ■ 

We now decompose the solution of f l9.20p into 

= p^ + Tq\ p-" EG, e A. (9.26) 

We have: 

Proposition 9.8 (i) is the solution of ^9.20^1 if and only if the pair {p'^,q'^) given by 
^9. 2(A) is the solution of the following variational formulation: 

Find [jf ,q'^) eV X A such that: 

a^{q\v) + {p^ + Tq\v)^= {T,v), \/v eV, (9.27) 
a^(p^,u') = 0, Vw e A (9.28) 

(ii) This variational formulation is well-posed and there exists a constant G, independent 
of T such that 

\{p\q^)\vxV<G\V\y,. (9.29) 

(Hi) When r — )■ 0, (p'^, g"^) — (p°, q^) and = p'^ + Tq'^ — j- p^ where {jp , q^) is the solution 
of the following variational formulation: 

Find (p°, g°) e V X A such that: 

a0(g°, v) + (p°, v)^ = (r, v), \fvE V, (9.30) 
a^{p^,w) = 0, yweA. (9.31) 



Proof: (i) We insert the decomposition fl9.26p into fl9.20p and get {p'^, q'^) E G x A such 
that (19.271) holds for any v E V. Now, using the characterization (19.251) of G, we get 
flOHD . Conversely, let (p^,g^) E V x A he a solution of flOTI) . flOHj) . Constructing 
according to (19.261) obviously leads to (I9.20p . 

(ii) We refer to [17|. 

(iii) is obvious. ■ 

Now, we examine what is the PDE problem solved by the variational formulation 
(I9.27p . (19. 28 p . Taking successively a smooth test function with compact support, followed 
by a smooth function with non-zero boundary values, we easily find that (p"^, q'^) solve the 
following PDE problem: 

e-^P^ + rgO - 4: ( ^-'^) = ^> [0, 1], (9.32) 



dx \ dx ^ 

(e-'^^m = g{0), (e-^^)(l) = ^(1), (9.33) 
ax dx 

q^iO) = 0, (9.34) 

p'^ is a constant. (9.35) 
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The boundary value problem for is overdetermined, having one Dirichlet and one 
Neumann boundary condition at a; = 0. However, the presence of the unknown constant 

introduces another degree of freedom which compensates for this over-determination. 

More precisely, let us consider the map J-" which to p'^ G M associates J^ip^) = {^~'^^~ 
g){0) G M, where q'^ is the unique solution of the mixed Dirichlet-Neumann boundary value 
problem: 

e-^p- + rgO " ^ fe"^^") = F, in [0, 1], (9.36) 



dx \ dx ^ 

g-(0) = 0, (e-'^^)(l)=^?(l), (9.37) 

This problem is uniquely solvable. Additionally, it is uniformly solvable when r — 
thanks to the Dirichlet boundary condition at x = which makes the elliptic operator 
— (e""^^) invertible with these boundary conditions. Finding a solution to (19. 32 p . 
(19.351) is equivalent to finding a root to the equation J-'(p^) = 0. But J^{j>^) is an affine 
function of with 

^(p)= J-(0)+^'(0)p, 
and J^'(O) = e~''^^^(0), with Q'^ the unique solution of the problem 

+ rQ^) - 4: (e-'^] = 0, in [0, 1], (9.38) 



dx \ dx ^ 

g^(0) = 0, (e-<^^)(l) = 0. (9.39) 

Then, the existence of p"^ is equivalent to -F'(0) being non-zero. By contradiction, sup- 
pose that -F'(0) = 0. This means that ^-(0) = 0. Then, denoting by = and 
differentiating f l9.38p with respect to x leads to 

tU^-^( e^-^(e-^f/")^ = 0, in [0, 1], (9.40) 
dx \ dx J 

t/-(0)=0, t/"(l) = 0, (9.41) 

which implies that = 0. It follows that Q'^ is a constant and by (19. 38 p . this constant 
must be Q'^ = — 1/t. But this is in contradiction with the first boundary condition (19. 39 p . 
Therefore, ^r(O) ^ and consequently J^'(O) 0. Finally, there exists a unique root p"^ 
of J^ip^) = given by p'^ = —yr^- Additionally, this root is uniformly bounded when 
r — J- 0, because both J-'(O) and J^'{0) remain bounded when r — )■ 0. Indeed, both elliptic 
problems associated to J^{0) and J^'(O) are uniquely solvable when r = 0, thanks to the 
Dirichlet boundary conditions and the Poincare inequality and J^'{0) cannot be zero also 
in that case. As a summary, we have proved by elementary means that problem (I9.32p - 
(I9.35P is uniquely solvable and that its inverse is uniformly bounded with respect to r 
when r — 0. 
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9.3 The multi-dimensional strongly anisotropic elliptic problem: 
the continuous case 



We recall some notations. We define h = B /\B\ where i? is a divergence free field {nabla- 
B = 0) over the domain Q. As a consequence, V ■ b = —b- V(ln|i?|). We define the 
boundaries r+, r_, Tq like in section 1631 The magnetic field lines are curves, solutions of 
the differential equation ^(s) = b{X{s)). For simplicity, we suppose that all magnetic 
field lines intersect the boundaries r+ and r_ at two points denoted by Xq and Xq 
respectively. In this case, we denote the magnetic field line starting from x^: X{xq;s), 
or simply B^~. If b is regular, by the Cauchy-Lipschitz theorem, for any x G fi, there 

exists a unique Xq such that x G B^-. We count the arclength s from its origin at Xq . 
Extensions of the present theory to more complex cases, where the magnetic field lines 
may intersect the boundary at one or zero points are possible but will be discarded here 
for simplicity. 

We introduce an 'arclength potential' by: 



0(3;) = - _Eiiiy)ds{y), (9.42) 



where Xq is such that x G B^- and Xq x denotes the arc along the magnetic field line 
issued from Xq and ending at x. Note that there is no condition on E for the existence of 
(j), because </> is not a 'true' potential, but only the primitive of the electric field along the 
direction of the magnetic field lines. The introduction of </> allows to remove the electric 
field by a simple transformation. Indeed, we remark that 

b ■ (Vn" - n^E) = e-^b ■ Wu\ = e V. (9.43) 

Then, we have the 



Proposition 9.9 With defined by { 9.4^ , problem AE is equivalent to the following 



modified anisotropic elliptic (MAE) problem: 



re-^u^ - V ■ (e'^ib ® b)Vu^) = rF, m n, (9.44) 
{b-u)e~%-Vu^ = Tg, onT+ur.. (9.45) 



We note the following lemma: 

Lemma 9.10 (i) Let u be a solution of the problem: 

V ■ {bu) = 0, m n, (9.46) 

{b-u)u = 0, on r+ ur_. (9.47) 

Then, u = identically on Q. 
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(a) Let f be defined on Q and g be defined on r+ U r_. Then, the problem 



V ■ (bu) = /, in Q, 

{b ■ v)u = g, on r+ U T^. 



(9.48) 
(9.49) 



has a solution provided the following solvability condition is satisfied: 



f 



9 



1 a:, 



.+ 




(9.50) 



- \B 



X, 







for all field lines, i.e. allx^ G r_. Under this condition, the system has a unique solution: 



where Xr. is such that x G . 

Proof: (i) Elementary computations show that b ■ V(^) = j^V ■ (bu). Then, ( I9.46P is 
equivalent to saying that is constant along each magnetic field line. But f l9.47p implies 
that u = 0. 

(ii) Eq. f l9.48p implies that ^ ■ V(^) = Integrating this equation along any 
magnetic field line leads to fl9.5ip . Identifying f l9.5ip at Xq with the boundary condition 
(I9.49P leads to the compatibility condition fl9.50p . If this relation is satisfied, then, formula 
(19.5 ip provides the unique solution to the problem. ■ 

With this lemma, we can now prove the 

Proposition 9.11 Let be the solution of the MAE problem. Then, when t 0, u"^ 
converges to where vP is constant along the magnetic field lines (i.e. b ■ Vm"^ = ^) O'nd 
is given by: 



where Xn is such that x G B- . 

Proof: We expand = vP + tv} + o(r). Inserting this expansion in the MAE problem, 
we get, at leading order: 




(9.51) 




(9.52) 



-V ■ {e^^{b®b)Wu^) = 0, in n, 
[b ■ v)e-'^b ■ Vu° = 0, on r+ u r_. 



(9.53) 
(9.54) 
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and at the next order: 

e-V- V- (e-'^(6®6)Vu^) = F, in (9.55) 
{h-u)e-%-S/u^ = g, onr+ur_. (9.56) 

From fl9.53p . fl9.54p and Lemma [9.101 (i). we deduce that h ■ VvP = 0. Then, from Lemma 
19.101 (ii). problem (I9.55p . (I9.56P is solvable if and only if the compatibility condition ( 19.50p 
is satisfied (with / replaced by F — e~'^u^). This leads to condition f l9.52p . ■ 

Back to the ri^ variable, we find that — )■ n° = u^{xQ)e~'^ with u'^{xq) given by f l9.52p . 

We now introduce a variational formulation. We denote by V the space V = {u E 
L^(f2) I b-Vu G L^(r2)} where L^{Q) is the space of square integrable functions on Q. We 
endow V with the norm \u\v = (|^^P + \b • Vnp)"*^/^ where \u\ = (J^ |u(x)p(ix)^/^ is the 
norm on L^(f2). The scalar products on V and L^(fi) are respectively denoted by {u,v)v 
and {u,v). We also define the following bilinear forms on V and L^(0, 1) respectively: 

a^{u,v)= / e~'^{b-Vu){b-Vv)dx, {u,v)4,= / e~'^uvdx. (9.57) 
Jn Jn 

Finally, we define the linear form F on ^ by: 

{r,v) = {F,v)+ gvdr{x) (9.58) 

Jr-ur+ 

Now, taking v E V, multiplying f l9.44p by v, integrating on x E Q and using the boundary 
condition (19.450 . we are led to the following 

Proposition 9.12 A variational formulation of problem MAE is as follows: 

Find eV such that: 

a^{u\v) + T{u\v)^ = T{V,v), MveV. (9.59) 



Like in the one- dimensional case, we note that the bilinear form a^{u, v) + t{u, v)rp is 
coercive on V. However, the best coercivity estimate is of the order of r and is given by 
(I9.2ip . This leads to an estimate of the solution m"^ in terms of the right-hand side F of the 
form fl9.22p . which deteriorates as r — )■ 0. Therefore, a standard finite element method 
based on ( 19.20p cannot be used when r ^ 1. We aim at finding a variational formulation 
which has bounded condition number when r — ?■ 0. For this purpose, like in the one- 
dimensional case we adopt the method of [I7j based on a micro-macro decomposition of 
the solution u'^. 

Definition 9.13 We define the spaces G and A as follows: 

G = {v EV\b-\/v = 0}, A = {v EV\v\r_ =0}. (9.60) 
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Proposition 9.14 (i) We have: 

V = G®A. (9.61) 

(ii) G is characterized by 

ueG ueV such that a^{u,v)=0, \/v e A. (9.62) 



Proof: (i) First, suppose that v & G (1 A. Then, v is constant along the magnetic field 
lines and its restriction on r_ is zero. By the assumption that all magnetic field lines 
contained in Q intersect r„ at one point, v is identically zero in which shows that 
G n A = ^. Let us now take any function v & V and construct a function p G G by the 
formula p{x) = v{xq), where is the foot of the magnetic field line passing at x (in other 
words, w e -^xq )• Then q = v — p obviously belongs to A. This shows that V = G + A. 

(ii): the proof of proposition 19. 71 can be reproduced in the mult i- dimensional case without 
any change. ■ 

Like in the one-dimensional case, we decompose the solution u'^ of (19.59P into 

=p^ + Tq\ G G, e A. (9.63) 

Then, proposition 19.81 applies without any change to the mult i- dimensional case. We state 
it in detail for the sake of completeness (the proof is similar and is omitted): 

Proposition 9.15 (i) is the solution of ^9.59i) if and only if the pair (p"^, q'^) given by 
^9. 63\) is the solution of the following variational formulation: 

Find {p'^ , q^) e V X A such that: 

a^iq\v) + {p^ + Tq\v)^= {T,v), E V, (9.64) 
a^{p^,w) = 0, \/weA. (9.65) 

(ii) This variational formulation is well-posed and there exists a constant G, independent 
of T such that 

\{p\qn\vxV<G\r\v'. (9.66) 

(Hi) When r — )■ 0, [p'^, q'^) — (p°, g°) and = p'^ + Tq'^ — )■ p° where (p°, q'^) is the solution 
of the following variational formulation: 

Find (p°, g°) G X A such that: 

a^iq", v) + (/, v)^ = (r, v), \fvE V, (9.67) 
a^{p^,w) = 0, WweA. (9.68) 
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Now, we derive the PDE solved by the variational formulation (19 .64^ . f l9.65p . Using the 
same method as in the one-dimensional case, we find that (p"^, q^) solve the following PDE 
problem: 



{b-v)e-^b-Vq^ = g, on r+ U r_, 

= on r„, 
-V ■ {e-^{b®b)Wp^) = 0, in 
{b ■ u)e-'t'b ■ S/p" = 0, on r+. 



in 



(9.69) 
(9.70) 
(9.71) 
(9.72) 
(9.73) 
(9.74) 



The boundary value problem for p"^ is equivalent to saying that p'^ is constant along the 
magnetic field lines (6 ■ Vp"^ = 0) and that p'^ is determined by its boundary value on 
r+. The boundary value problem for q'^ is overdetermined, having one Dirichlet and one 
Neumann boundary condition on r_. However, the presence of the unknown constant 
p'^\r+ introduces just the right number of degrees of freedom to compensate for this over- 
determination. 

We can also see intuitively that this formulation has a condition number which is inde- 
pendent of r. Indeed, for given p"^, constant along the magnetic field lines and determined 
by its boundary value on r_|_, we can find q'^ by solving the mixed Dirichlet-Neumann 
boundary value problem: 



The operator q — )■ re'^q — V ■ {e~'^{b ® b)'Wq) with Dirichlet boundary conditions 
on r_ and Neumann boundary conditions on r_|_ is invertible whatever the value of r 
is. Therefore, the inversion of (19. 75p - ( 191771) provides a solution q^ as a function of p"^, 
whatever the value of r. The additional boundary condition {b ■ i')e~'^b ■ Vg"^ = g on 
r_ gives another condition which allows us to compute p"^. The solvability of this last 
condition is also indpendent of the value of r as was explicitly seen in the one-dimensional 
case. 

9.4 The discrete anisotropic elliptic problem: back to the one- 
dimensional problem 

The discretization of the Euler-Lorentz model directly leads to an anisotropic elliptic 
problem in discrete form. Therefore, we are not free of choosing the numerical method 
for dicretizing it. Rather, the method is imposed by that of the discretization of the 
background Euler-Lorentz model. Despite this fact, we will see that we can apply the 
same ideas than those developed in the last two sections for the continuous problem, to 



e-'^(p" + rg") - V ■ (6-^(6 ® 6)Vg") = F, 
[b ■ iy)e~^b ■ Vg" = g, on r+, 
g"^ = on r_. 



in 



(9.75) 
(9.76) 
(9.77) 
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these discrete formulations. Like for the continuous problem, we will start to develop the 
ideas in a simple one- dimensional framework, before applying them in full generality to 
the three (or more)-dimensional framework. 

The discretization of the Euler-Lorentz model by the FDAP-2 scheme (a similar study 
could be conducted for the FDAP-1 scheme) leads to a discretization of the form: 

[nk+2 -riK) - tik+iEk+i 



2 

{riK - nK-2) - uk-iEk-i^ I + tuk = tEk, (9.78) 

where again, we restrict ourselves to linear pressure-density relationships. We also have 
made 5 = 1 and collected all known terms at the right-hand side into the generic term 
Fk- K is now a one-dimensional index ranging in Z. The discrete electric field is also 
one- dimensional and denoted by Ek- Eq. fl9.78p is a generic one- dimensional model for 
( 18391) . 

Now, we introduce the boundary conditions which are an important aspect of this 
discussion. We suppose that the unknown rtx is defined for K ranging over a finite 
interval K e [-M, M] := {-M, -M + 1,...,M-1, M}. For the sake of the forthcoming 
developments, we write ( I9.78p in the form of a mixed problem and we highlight the 
dependence upon r: 

(7x+i - Yk^,) + rnl- = tEk, K G [-M, M], (9.79) 

7^ = ^K+i-^^_i] -^^^i^, Ke[-M + l,M-l]. (9.80) 

To complete this formulation, we need to impose boundary conditions in the last row of 
cells {K = M or K = —M) and in an additional row of fictitious cells {K = M + 1 or 
K = —M — 1). We impose the following boundary conditions (which are the discrete 
counterpart of f l6.15p ): 

7m = 7m+i = t-Qm, Y^m = Ym~i = (9.81) 

where g±M is supposed known. 

Since our solution method for strongly anisotropic elliptic problems relies on the intro- 
duction of an appropriate variational formulation, we first introduce the discrete counter- 
part of the original variational formulation fl9.20p . Multiplying fl9.79p by a test sequence 
{vk}k£[-m,m]^ performing a discrete integration-by-parts and using fl9.8Qp together with 
the boundary conditions ( 19.8ip . we get the following 

Proposition 9.16 The solution n]^ of problem \9. 79\ ), l{9.8(]\) with boundary conditions 
Ii9.81\) satisfies the following discrete variational formulation: 

Find n'^ eV such that: 

aE{n\v) + T{n\v) = T{T,v), Vt; G V, (9.82) 
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where V = i^{[—M,M]), the space of square integrable sequences on [—M,M], and 

aE{u, v) = ^ — I — [uk+i - uk~i] - ukEk \ {vk+i - vr-i) , (9.83) 

k=-M+l ^ ^ 

M 

(u,v) = ^ UkVk, (9.84) 

k=-M 

(r, v) = {F, v) + — [gM{vM^i + vm) - g^M{v-M+i + V-m)] ■ (9.85) 



We note that the form qe is not symmetric but its leading order part (corresponding 
to the discretization of the second order derivative) is. We also note that the for aE{u, u) + 
t{u, u) is coercive on V . However, we have no better estimate than 

aE{u,u) + t{u,u) > r\u\y, (9.86) 

(with l^ly = {u,u)) which implies that there is no better estimate for the solution n'^ of 
( I9.82P than \n'^\v < ''"~^|r|v. The condition number of formulation f l9.82p tends to oo as 
r — 7- and it cannot be used for solving for n'^ when r is very small. We are going to 
remedy to this problem by introducing the discrete analog of formulation fl9.27p . f l9.28p . 

Definition 9.17 We define the spaces Ge and A as follows: 

GE = {veV \^ K+i - VK-i] - vkEk = 0, VA' G [-M + 1, M - 1]}, (9.87) 
A = {veV\ v_M = v_M+i = 0}. (9.88) 



Proposition 9.18 (i) We have: 

V = Ge®A. (9.89) 

(ii) Ge is characterized by 

u^Ge ueV such that aE{u,v) = 0, \/v e A. (9.90) 



Proof: (i) Let v G Ge A. Then, vk is the solution of a two-stages linear recursion 
where two successive steps are zero. This implies that v is identically zero, proving that 
Ge A = 0. Now, let v & V. Then, there exists a unique p & Ge which satisfies 
P-M = V-M, P-M+i = v_M+i- Indeed, p is the solution of two-stages linear recursion 
and has specified values at two successive points. Then, px is uniquely defined for all 
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K's. Then, defining Qk = vk — Pk, we obviously have g G A, showing that V = Ge + A 
and ending the proof of point (i). 

(ii): that u & G =^ aEiu.v) = 0, Vt" G A is obvious (it is even true for all v G V). 
Conversely, suppose that u E V such that aE{u,v) = 0, Vf G A. Let ipx for K G 
[— M + 1, M — 1] be arbitrary. There exists v E A such that ipx = vk+i — vk-i- Indeed, 
vk is defined by a two-stages recursion with initial conditions V-m = ^-m+i = 0. Then, 
such a. V in A exists and is unique. Replacing vk+i — vk-i by ipx in asiu, v) = 0, we are 
led to the following relation: 

— f — [^^fc+i - Uk~i] - UkEk j = 0, V{V'i^}Ke[-M+i,M-i]- 

This clearly implies that u & Ge and ends the proof of the proposition. ■ 
We now decompose the solution of fl9.82p into 

7f =p^ + Tq\ G Ge, G A. (9.91) 

We have: 

Proposition 9.19 (i) rf is the solution of Ii9.82\) if and only if the pair (p^, g^) given by 
l[9.91\) is the solution of the following variational formulation: 

Find {p'^ , q'^) e V X A such that: 

aE{q\v) + ip^ + Tq\v) = {T,v), e V, (9.92) 
aE{p^,w) = 0, ^weA. (9.93) 

(ii) This variational formulation is well-posed and there exists a constant G, independent 
of T and h such that 

\{p\qn\vxV<G\r\v^. (9.94) 

(Hi) When r — )■ 0, (p^, g^) — (p*^, g°) and = p'^ + Tq'^ — )■ p^ where (p°, g°) is the solution 
of the following variational formulation: 

Find [jp, q^) x A such that: 

aE{q\v) + {p\v) = {T,v), \/v G V, (9.95) 
a^{p°,w) = 0, yweA. (9.96) 



Proof: (i) We insert the decomposition ( I9.9ip into f l9.82p and get (p^, g^) & Ge x A such 
that fl9.92p holds for any v ^ V. Now, using the characterization fl9.90p of Ge, we get 
flCTj) . Conversely, let {p^,q^) G \^ x A be a solution of flCTj) . f lCTj) . Constructing ra^ 
according to fl9.9ip obviously leads to (19. 82 p . 
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(ii) We refer to [H]. 

(iii) is obvious. 



Another, more elementary view of this proposition is as follows. A given E Ge depends 
on two arbitrary quantities: pL^/ and p'^_j^.j j^i- For a given G Ge-, we can solve for 
satisfying f l9.92p . The form ao(M,w) is coercive on V. Indeed, if aoiu.u) = 0, then 
u E Go n A and with ( I9.89p . is such that u = 0. Additinally, the coercivity constant can 
be proven uniform with respect to h. Therefore, using test sequences v E A, f l9.92p can 
be solved uniquely for G A. Then, taking a sequence v E V such that vk = 0, for all 
K G [— M + 2, M] leads to two additional linear relations which allow to determine p'Lm 
and p'Lm+1 uniquely. This procedure allows us to determine {p'^, q^) uniquely in a uniform 
way with respect to both r and h. 

We now turn ourselves to the multi-dimensional case. 



9.5 The discrete strongly anisotropic elliptic problem: the multi- 
dimensional problem 

We now assume that the multi-index K belongs to a box /C = Mj, Mj]. The 

discretization of the Euler-Lorentz model by the FDAP-2 scheme (a similar study could 
be conducted for the FDAP-1 scheme) leads to a discretization of the form: 



i,j=l ^ ^ ' 



—bj\K-ejbi\K-ej ^~2~ ~ ^\K-ej-ei] — n\K-ejEi\K-e-^^ 

+Tn\K = tF\k, (9.97) 

where we have introduced a linear pressure-density relationships. We also have made 
6 = 1 and collected all known terms at the right-hand side into the generic term F\k- 

Now, we write f l9.97p in the form of a mixed problem and we highlight the dependence 
upon r: 

i=i 

Y\k = J2^^\K[^ in^\K+e. ' n^lx-e,) ' U^WE^W J , V/C G /Ci,t, (9.99) 

i=l ^ ^ 

where /Cint denotes the internal cells /Cint = nLi[~"^« + l,iVfj — 1]. To complete this 
formulation, we need to impose boundary conditions in the last row of cells {K E ICj = 
nLi[-^i' ^i] \ U!i=i[-^i + 1, Mi - 1]) and in an additional row of fictitious cells {K G 
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ICf = YlLi[-Mi - l,Mi + 1] \ Yll^-^[-Mi,Mi]). We impose the following boundary 
conditions (which are the discrete counterpart of (16.151) ): 



^jkIk Ik = TdK, G /C/ U /Ci? and jx such that uk = ±e 



(9.100) 



where qk for K ^ ICj U JCp is supposed known. The condition on jx is that the direction 
Jk is one of the directions to which the boundary of /C at cell K is normal. 

Multiplying (I9.98P by a test sequence {vk}kgic, performing a discrete integration-by- 
parts and using (I9.99P together with the boundary conditions (I9.100p . we get the following 

Proposition 9.20 The solution n]^ of problem Ii9. 98\) . i\9. 99\) with boundary conditions 
l[9.100\) satisfies the following discrete variational formulation: 



Find n'^ eV such that: 
where V = (?{JC), the space of square integrable sequences on fC, and 



(9.101) 



aE[u,v 



E 



E"^ 



•i K 



K 




{v\K+e, - v\K-e,) 



E 



-1 



'3K 



is defined by uk 



(9.102) 
(9.103) 

(9.104) 



{T,v) = iF,v) + 

where in the last term defining {T,v), 

We again note that aE{u,u) + t{u,u) is coercive on V. However, we have no better 
estimate than fl9.86p . which implies that the condition number of formulation fl9.10ip 
tends to oo as r — 0. Therefore, this formulation cannot be used for solving for n'^ when 
r is very small. We then introduce a new formulation in the spirit of what has been done 
in the one- dimensional case. To this aim, we assume that there is a space direction io 
such that big\K for all cells K E K. This condition can be removed at the expense of 
some additional work which will be detailed in future work. We suppose that = 3 for 
simplicity. We also define the set Kg by: 

3 2 

1Cg = \{[-M,,M,] \ J][-M, + 1,M,-1] X [-M3 + 2,M3]. (9.105) 



i=l 



1=1 



We note that in the 3rd direction (corresponding to Iq in the general case), the boundary 
cells are shifted and two layers of boundary cells are defined on the left-hand side and no 
boundary cells on the right-hand side. The complement /C \ Kg = Y\d=i[~-^i + 1, — 
1] X [-Mg 2, Ms] is denoted by Ka- 



87 



Definition 9.21 We define the spaces Ge and A as follows: 

3 

^bi\K {i2hi)~^ iu\K+ei - u\K-e,) - M|A'-Ei|x) = 



Ge= <veV 



i=l 



A = {v \v\k = 0,yK e /Cg}. 



\/K eJCm}, (9.106) 
(9.107) 



Proposition 9.22 (i) We have: 

V = Ge®A. (9.108) 

(a) Ge is characterized by 

ueGe ueV such that aE{u,v) = 0, \fv e A. (9.109) 



Proof: (i) We note that, for the elements p of Ge, the components p\k for K G Kg are 
free. Similarly, for any element q & A, the components qlx G ICa are free. Now, Let 
V G n A. Because v E A, v\k is identically zero on JCg- Then it is easy to see that the 
condition that v E Ge, owing to the fact that 63 |a' 7^ 0, allows to determine v recursively 
in the layers corresponding to K3 = Constant, starting from the layer K3 = — M3 + 2 
up to the layer K3 = M3. And this recursive calculation shows that v\k is identically 
zero on /C. This shows that Ge A = 0. Now, we consider v E V. Using the same 
recursive procedure as above, we can define an element p G Ge such that p\k = v\k for 
all K G JCg- Once p is found, q = v — p is defined and obviously belongs to A, showing 
that V = Ge + A. 

(ii): that u E G =^ aE{u,v) = 0, Wv E A is obvious (it is even true for all v G V). 
Conversely, suppose that u eV such that aE{u,v) = 0, Vf G A. Let ipK for K G /Cint be 
arbitrary. There exists v E A such that 



(9.110) 



Indeed, vk is defined by the same recursive procedure as above, starting from the zero 
values of v\k for K E Kg- Additionally, the so-constructed v is the unique v E A 
satisfying property (19.1101) . Inserting (19.1101) into the expression (I9.102p of a^; and using 
that aE{u,v) = 0, we are led to the following relation: 



E 



KefC: 



int L 



i=l 



\K 



[U\K+e 



U K- 



K 



K 



0. 
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This clearly implies that u G Ge and ends the proof of the proposition. ■ 
We now decompose the solution n'^ of flQ.lOip into 

n^ = p^ + Tq\ p^eGs, g"eA (9.111) 
We can copy proposition 19. 191 'mutatis mutandis': 

Proposition 9.23 (i) is the solution of Ii9.101\) if and only if the pair {p'^,q'^) given 
by Ii9.111\) is the solution of the following variational formulation: 

Find {p'^ ,q'^) e V X A such that: 

aEiq\v) + ip^ + Tq\v) = {T,v), Wv e V, (9.112) 

aE{p^,w) = 0, ^weA. (9.113) 

(a) This variational formulation is well-posed and there exists a constant G, independent 
of T and h such that 

\{p\qn\vxV<G\r\v:. (9.114) 

(Hi) When t ^ 0, (p^, g^) — (p°, g°) and = p'^ + rq'^ p^ where (p°, g°) is the solution 
of the following variational formulation: 

Find (p°, q^) e V X A such that: 

aE{q\ v) + (/, v) = (r, v), yv G V, (9.115) 

a^{p°,w) = 0, VweA (9.116) 



This variational formulation has the following interpretation: First, take an arbitrary 
G Ge- Then, depends on independent components p\k for K belonging to JCq- 
Using the fact that a^; is coercive on A (the proof is similar as in the one-dimensional 
case), we find a unique q"^ E A such that fl9.112p holds for any v G A. But q'^ depends 
on the chosen element of G. Then, taking v such that t'l;^ 7^ if and only ii K E JCq 
gives as many independent relations as needed to fully determine p'^. 

Variational formulation f l9.112p . fl9.113p leads to a solution methods for the discrete 
anisotropic elliptic problem with a uniform condition number with respect to r when 
r ^ 1. This method does not require the numerical determination of the magnetic field 
lines nor any integration along these lines. 

9.6 Strongly anisotropic elliptic problems: conclusion 

The resolution of the Euler-Lorentz model when r ^ 1 (the limit r — > corresponding 
to the so-called drift-fluid limit) leads to a strongly anisotropic discrete elliptic problem. 
We have seen that this problem degenerates when r — )■ and leads to an ill-conditionned 
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numerical resolution. We have proposed a new variational method. This method has 
a condition number which is independent of r when r <C 1 and is therefore efficient 
independently of the value of r. Additionally, it provides the correct solution to the 
limit problem when r — 0. The knowledge of the magnetic field lines is not needed 
and no integration along these field lines need to be performed. Therefore, this method 
is particularly suitable to a context where the magnetic field is susceptible to vary with 
time. 
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Conclusion 



In these notes, we have described how to construct Asymptotic-Preserving schemes for 
plasma fluid models in a variety of situations. We have first considered the quasi-neutral 
limit and applied the methodogology to the Euler-Poisson and to the Euler-Maxwell 
problems. In a second part, we have focused on the Euler-Lorentz model in the drift-fluid 
limit which arises when the magnetic field is large and simultaneously the Mach number is 
small. In all cases, the same methodology has been applied. First, we find a reformulation 
of the original problem in such a way that it directly appears as a perturbation of the 
limit problem. Then, we focus on the time discretization by considering time semi-discrete 
schemes and we determine which terms must be evaluated implicitly in order to make the 
scheme AP. Once a proper time discretization has been found, we apply it to a fully- 
discrete version of the scheme. In this last step, issues like conservativity or numerical 
viscosity can be brought into the framework of AP schemes. 

We have focused on a presentation of the methodologies. We refer to the bibliography 
given in the introduction for applications to practical cases and performance tests. 

The AP methodology can be applied to a large variety of situations. Of particular 
interest are cases where where several limits must be taken independently. An important 
issue for instance in two-fluid models is to treat simultaneously the smallness of the Debye 
length (quasi-neutral limit) and the smallness of the electron mass (Low Mach-number 
limit in the electron fluid equations). Another issue is the design of schemes for the Euler- 
Lorentz model which are AP when both the Debye length and the cyclotron period may 
tend to zero independently. 

Other open problems concern the stability analysis of the schemes in the nonlinear 
settings, as well as the obtention of rigorous error estimates. 
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